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The  National  Bureau  of  Standards  has  been  engaged  in  the  heat 
transfer  analysis  of  underground  installations  for  the  past  several 
years.  This  report  covers  a part  of  the  Bureau’s  activities  related  to 
the  computer  simulation  of  the  thermal  environment  for  prototype  shelters* 

The  computer  was  used  basically  to  simulate  energy  balance  in 
the  shelter  living  space  and  to  analyze  heat  conduction  from  the  shelter 
walls  (including  ceiling  and  floor)  to  the  surrounding  earth. 

For  the  heat  conduction  analysis,  finite  difference  techniques 
were  employed;  using  a three  dimensional  model  in  some  cases  and  a one 
dimensional  model  for  the  remainder* 

Digital  computer  programs  were  developed  and  applied  to  seven  different 
prototype  shelters  for  which  temperature  and  humidity  records  with  simulated 
occupants  were  available  as  a result  of  studies  by  the  National  Bureau  of 
Standards  and  by  the  University  of  Florida.  In  the  seven  shelters  used 
for  the  investigation,  twelve  different  operating  conditions  were  analyzed. 

Of  these  twelve  conditions,  ten  were  under  summer  operation  and  two  under 
moderate  winter  conditions. 

Generally  the  agreement  between  the  computed  and  observed  thermal 
environment  on  these  prototype  shelters  was  surprisingly  good,  in  spite  of 
the  fact  that  numerous  simplifications  were  involved  in  describing  the  complex 
shelter  heat  transfer  system  for  computer  analysis.  Two  inherent  uncertainties 
exist,  which  influence  the  final  reliability  of  the  calculations.  The  first 
involves  the  description  of  the  actual  complex  system  by  mathematical  language 
(or  operational  uncertainty).  The  second  is  related  to  the  accuracy  of  input 
data  used  for  the  calculations,  (or  data  uncertainty).  Often,  these  are 
interrelated. 
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In  this  study,  the  following  four  different  computer  models  were 
studied  and  two  of  them  were  extensively  utilized  for  the  comparison 
of  calculated  thermal  environment  with  the  observed  data  in  the  proto- 
type shelters. 


M-(l) 

Three-dimensional  rectangular  model  with  composite  walls 
and  with  separate  initial  temperature  patterns  normal  to 
the  six  bounding  surfaces. 

M-(2) 

Three-dimensional  rectangular  model  with  homogeneous  heat 
conduction  medium  having  separate  initial  temperature 
patterns  normal  to  the  six  boundary  surfaces. 

M-0)  One-dimensional  compound  model  for  six  composite  wall 
systems. 

M-  (4)  One-dimensional  compound  heat  conduction  model,  same  as 
M-0)  except  that  the  roof  region  was  assumed  adiabatic. 

To  simulate  the  initial  earth  temperature  distribution,  the  follow- 
ing three  modes  were  employed: 


i-d) 

Earth  temperature  gradients  normal  to  the  six  bounding 
surfaces . 

I-(2) 

Earth  temperature  gradient  normal  only  to  the  ground 
surface . 

l-(3) 

Initial  earth  temperature  constant  around  the  shelter. 

One  of  the  factors  not  well  established  for  calculating  the  shelter 
heat  transfer  is  the  heat  exchange  between  the  shelter  air  and  the  inner 
surfaces,  between  the  occupants  and  the  surfaces,  and  among  the  surfaces 
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The  analysis  of  simultaneous  exchange  for  radiative  and  convective 
energy  among  occupants,  air  and  inner  surfaces  of  a shelter  is 
very  complex,  and  it  requires  the  solution  of  a set  of  integral 
equations  which  are  difficult  to  solve  for  even  very  simple  geo- 
metries. Therefore  conventional  combined  heat  transfer  coefficients 
for  radiation  and  convection  were  used  in  the  analysis.  Several 
numerical  values  and  combination  of  these  combined  coefficients 
were  assigned  to  the  six  interior  surfaces  to  study  the  overall 
effect  on  shelter  thermal  environment. 

Findings 

1)  A soil  analysis  of  the  earth  around  most  of  the  prototype  shelters 
indicated  that  the  thermal  diffusivity  and  thermal  conductivity 

2 p 

were  in  the  neighborhood  of  0.02  ft  /hr  and  0.75  Btu/hr,  (ft)  , 
°F/ft  respectively.  These  values,  in  turn,  seem  to  result  in 
a good  agreement  between  the  calculated  and  the  observed  earth 
temperature  change  surrounding  prototype  shelters. 

2)  The  following  combined  heat  transfer  coefficients  at  the  shelter 
inner  surfaces  produced  satisfactory  simulation  of  the  shelter 
summer  environment  for  most  of  the  prototype  shelters. 

1.0  Btu/hr.  (ft1 2),  (°F)  for  vertical  walls 
1.5  Btu/hr.  (ft2),  (°F)  for  the  ceiling 
0.5  Btu/hr.  (ft2),  (°F)  for  the  floor 
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Although  there  may  be  some  other  values  and  other  combinations  of 
these  values  that  might  have  resulted  in  a slightly  better  simulation 
than  those  ussd  in  this  analysis,  these  three  values  can  be  considered 
representative  design  heat  transfer  coefficients  in  the  underground 
cavities . 

3)  For  larger  shelters,  the  one -dimensional  and  compound  model  (M-Q)) 
will  probably  be  adequate  for  calculating  the  shelter  thermal 
environment.  The  complicated  three-dimensional  model,  therefore, 
may  not  be  required  for  the  calculation  simulating  the  14-day 
occupancy  of  many  large  community  shelters.  For  small  shelters 
(such  as  family  shelters  similar  to  the  NBS  shelter) , however,  it 
is  recommended  that  the  three-dimensional  model  be  used  for  the 
accurate  calculation. 
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NOMENCLATURE 


Unless  otherwise  defined  in  the  text,  the  symbols  used  in  this 
summarized,  as  follows: 

report  are 

Symbols 

Dimensions 

a 

shelter  dimension  (half  of  the  inside  length) 

ft 

b 

shelter  dimension  (half  of  the  inside  width) 

ft 

C 

P 

specific  heat  of  moist  air 

Btu/lb, °F 

c 

shelter  height 

ft 

d 

shelter  depth,  distance  between  the  earth  surface 
and  the  ceiling  of  the  shelter 

ft 

G 

flow  rate  of  ventilation  air 

CFM 

h 

g 

ground  surface  heat  transfer  coefficient 

Btu/hr,  ft2,°F 

h or  h 
K 

k 

inner  surface  heat  transfer  coefficient  of  K th 

Btu/hr,  ft2,  °F 

exposure 

thermal  conductivity  of  solid 

Btu/hr,  ft,  °F 

le 

* 

Lewis  Relation  - _JL_ 

crC 

P 

dimensionless 

Pvs 

saturated  vapor  pressure  of  water  at  temp,  t 

inches  Hg 

Pv 

vapor  pressure  of  water  in  the  air  at  dew  point 
temperature 

inches  Hg 

PB 

barometric  pressure 

inches,  Hg 

QVS 

sensible  heat  released  by  the  ventilation  air 

Btu/hr 

qvl 

latent  heat  released  by  the  ventilation  air 

Btu/hr 

qgs 

sensible  heat  generated  in  the  shelter  by 
simulated  occupants 

Btu/hr 

Qgl 

latent  heat  generated  in  the  shelter  by 
simulated  occupants 

Btu/hr 

qwsk 

sensible  heat  released  by  the  shelter  inner  surface 
of  K th  exposure 

Btu/hr 

qwlk 

latent  heat  released  by  the  shelter  inner  surface 
of  K th  exposure 

Btu/hr 

qms 

sensible  heat  generated  in  the  shelter  by  things 
other  than  simulated  occupants 

Btu/hr 
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Symbols 


Dimensions 


%L 


latent  heat  generated  in  the  shelter  by  things  Btu/hr 

other  than  simulated  occupants 


QSUN 

SK 


t 


W 

a 


solar  radiation  intensity  at  the  earth's  surface 
inner  surface  area  of  K th  exposure 
temperature 
humidity  ratio 


Btu/hr, ft2 

ft2 

°F 

lb  of  water  vapor/ 
lb  dry  air 


W 


s 


humidity  ratio  of  air  saturated  by  water  vapor  lb  of  water  vapor/ 

lb  of  dry  air 


(K  =s  1,6)  coordinate  system  used  for  shelter  heat  transfer 
q (K  = 1,6)  coordinates  of  the  system  boundaries 


K 


a 

T] 

0 

X 

a 

P 

AXK 

A9 

2 


thermal  diffusivity  of  earth  in  the  region  sur- 
rounding K th  exposure 

albedo  of  earth  surface 

time  coordinate 

latent  heat  of  vaporization  of  water 

water  vapor  transfer  coefficient  (lb/ (hr) (ft2 ) (lb/lb 

density  of  moist  air 

finite  difference  length  along  X^ 

finite  difference  time 

summation  symbol 


ft 

ft 

ft2 /hr 

dimensionless 

hr 

Btu/ lb 
dry  air) 
lb/ ft3 
ft 
hr 


Subscripts 

Unless  otherwise  stated,  the  following  rules  of  subscripting 
will  apply  to  all  of  the  variables, 
a shelter  space  properties 

outdoor  air  properties 


o 


V 


Subscripts-  - continued 


v 


ventilation  air  properties 


c 


concrete  property 


g 


ground  surface  properties 


00 


deep  underground 


K innersurface  exposure  index 

K = 1 = North 

2 = South 

3 = East 

4 = West 

5 **  Floor 

6 = Roof 

In  some  cases,  the  subscript  W is  used  to  denote  the  wall  properties  instead 
of  K being  1,  2,  3,  and  4. 

The  subscripts  R and  F are  employed  in  the  same  manner,  denoting,  respec- 
tively, the  properties  pertaining  to  roof  and  floor  regions. 

S sensible  heat  property 

s saturated  air  property 

L latent  heat  property 


Operational  Symbols 
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Fig . 
Fig. 

Fig. 

Fig. 

Fig. 

Fig. 

Fig. 

Fig. 

Fig. 
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COMPARISON  OF  DIGITAL  COMPUTER  SIMULATIONS  OF  THERMAL  ENVIRONMENT 
IN  OCCUPIED  UNDERGROUND  PROTECTIVE  STRUCTURES  WITH  OBSERVED  CONDITIONS 


By 

T.  Kusuda  and  P.  R.  Achenbach 
1.  INTRODUCTION 

This  report  compares  digital  computer  calculations  of  thermal  environ- 
ment for  underground  protective  shelters  with  observed  conditions  of  tem- 
perature and  humidity.  Several  digital  computer  programs  have  been 
developed  by  the  National  Bureau  of  Standards  for  the  purpose  of  simulating 
the  heat  transfer  of  underground  structures . These  computer  programs  were 
applied  to  7 shelters,  whose  thermal  environment  under  simulated  conditions 
of  occupation  had  been  observed  experimentally.  Since  the  thermal  environ- 
ment in  underground  protective  structures  may  become  extremely  unfavorable, 
particularly  during  the  summer  occupancy  period,  for  large  areas  of  the 
United  States,  the  majority  of  prototype  shelters  mentioned  herein  were 
tested  under  summer  climatic  conditions . Of  the  7 shelters  whose  thermal 
environments  were  calculated  and  compared  with  experimental  observations, 

12  different  test  conditions  were  included,  2 of  which  were  under  moderate 
winter  conditions . 

Analytical  and  experimental  studies  of  various  shelters  have  shown 
that  the  temperature  and  humidity  within  the  occupied  underground  shelter 
depend  on  many  parameters,  which  may  be  classified  as  follows: 

1 . Structural  characteristics  » 

a.  Size  and  shape. 

b.  Physical  and  thermal  properties  of  construction  material. 
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2 . Site  characteristics. 

a.  Physical  and  thermal  properties  of  earth  surrounding  the 
she  Iter . 

c.  Thickness  of  earth  cover. 

d.  Type  of  earth  surface  and  landscape. 

e.  Neighboring  buildings  and  installations. 

3 . Climatic  factors. 

a.  Earth  temperature. 

b.  Psychrometric  condition  of  outdoor  air. 

c.  Solar  radiation. 

d.  Precipitation. 

4 . Operational  characteristics . 

a.  Ventilation  rate. 

b.  Psychrometric  condition  of  ventilation  air 

c.  Density  of  Occupancy. 

d.  Activity  of  Occupancy. 

e.  Heat  and  moisture  release  by  equipment  in  shelter. 

f.  Emergency  condition  such  as  sealed  up  or  surface  fire  conditions. 
Testing  of  underground  structures  to  cover  even  a small  portion  of  all 

of  the  possible  combinations  of  the  above  parameters  is  a formidable  and 
expensive  task.  However,  the  number  of  tests  could  be  drastically  decreased, 
and  the  efficiency  of  testing  improved,  if  the  effect  of  various  parameters 
in  the  thermal  environment  of  a shelter  could  be  predicted  by  computation. 

The  mathematical  formulation  of  such  a computation  should  take  into  account 
a majority  of  the  important  parameters  so  the  sensitivity  of  the  overall 
thermal  environment  to  the  several  parameters  could  be  studied  individually 
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or  simultaneously  with  others.  The  mathematical  procedures  should  be 
simple  enough  so  the  computation  time  (or  computer  cost)  would  be  reason- 
able. Finally,  and  most  important,  the  computed  results  should  be 
re  liable . 

Some  previous  computations  of  the  thermal  environment  in  shelters 
have  been  reported  which  take  into  account  simulated  human  metabolism, 
ventilation  effects,  and  heat  transfer  in  the  earth  [1,2,3],  but  there 
have  been  very  few  actual  comparisons  between  calculations  and  observa- 
tions for  a given  system  over  a substantial  period  of  time.  It  is  for- 
tunate that  the  observed  results  of  six  of  the  seven  shelters  covered  by 
this  study  were  so  well  documented  by  reports  of  the  University  of 
Forida  [4] , thus  making  possible  comprehensive  comparisons  between  the 

calculated  thermal  environment  and  the  observed  results. 

2.  BASIC  HEAT  TRANSFER  RELATIONS 


Since  the  details  of  the  numerical  technique  employed  for  the  heat 
transfer  of  underground  protective  structures  have  been  reported  pre- 
viously [ 1]  , only  basic  mathematical  formulations  employed  for  all  of 
our  computer  programs  are  given  here: 

2.1.  Shelter  air  heat  balance. 
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2.2.  Shelter  Inner  Surfaces 
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2.3.  Concrete  (or  inner  wall)  heat  conduction 
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2.4.  Boundary  between  the  concrete  (or  inner  wall),  and 
earth. 
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2.5.  Earth  heat  conduction. 
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2.6.  Earth  boundary  conditions. 
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This  solar  heat  model  is  essentially  the  same  as  the  sol-air  tempera- 
ture concept  and  ignores  the  direct  radiation  heat  exchange  between  earth 
surface  and  sky,  but  it  does  include  the  effective  radiation  heat  exchange 
between  the  earth  surface  and  ambient  air  by  adjusting  the  value  of  N . 

g 

According  to  the  last  equation,  however,  the  earth  surface  temperature 
never  becomes  lower  than  the  air  temperature,  which  is  not  always  the  case. 

2.7.  Psychr ometr ic  calculations. 

Taking  advantage  of  the  large  memory  of  the  high  speed  digital  computer 
of  the  National  Bureau  of  Standards,  all  the  psychrometric  ca i -ulations  were 
performed  using  the  thermodynamic  properties  of  moist  air  published  by  Goff 
and  Gratch . 

The  thermodynamic  properties  of  dry  air  and  those  of  saturated  air  at 

one  standard  atmospheric  pressure,  such  as  the  following,  were  tabulated  as 

temperature  functions  by  Goff  and  Gratch  [ 5]  . 

Wg  = humidity  ratio  of  the  saturated  moist  air  (lb/lb  of  dry  air). 

h = enthalpy  of  dry  air  (Btu/lb  of  dry  air). 

hg  ~ enthalPy  °f  the  saturated  moist  air  (Btu/lb  of  dry  air)c 

hw  = enthalpy  of  the  water  (Btu/lb  of  water). 

Va  = volume  of  the  dry  air  (cu  ft/lb  of  dry  air). 

Vg  = volume  of  the  saturated  moist  air  (cu  ft/lb  of  dry  air). 

fg  = factors  related  to  the  relative  humidity  and  degree 
(dimensionless)  of  saturation. 

These  properties  were  read  into  the  computer  for  the  temperature  range 
from  30  °F  to  120  °F  at  every  one  degree  increment,  except  that  the  humidity 
ratio,  Ws,  was  programmed  for  the  temperature  range  from  -20  °F  to  120  °F 
at  every  one  degree  increment . 
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The  following  psychrometr ic  symbols  and  formulas  are  used  to 
derive  the  desired  properties  from  given  sets  of  properties,  such  as  dry 
and  wet-bulb  temperatures: 


Pv  = 

: partial  water  vapor  pressure  in  moist  air  (in.  Hg). 

PB  ‘ 

: barometric  pressure  (in.  Hg) . 

T) 
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II 

: partial  water  vapor  pressure  in  saturated  air  (in.  Hg . ) 

9 = 

: relative  humidity,  as  a fraction 

M* 

: degree  of  saturation. 

W 

: humidity  ratio  of  moist  air  (lb/lb  of  dry  air). 

V 

: volume  of  moist  air  (cu  ft/lb  of  dry  air). 

Ws*  = 

: Wg  evaluated  at  the  thermodynamic  wet-bulb  temperature 
(lb/lb  of  dry  air). 
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: h^  evaluated  at  the  thermodynamic  wet -bulb  temperature 
(Btu/lb  of  water). 
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h = p, (hg  - ha)  + ha  2.7-6 

h + (Ws*  - W)^*  = hg*  2.7-7 

The  last  formula  represents  the  thermodynamic  wet-bulb  temperature  relation 
for  given  h and  W of  moist  air. 

The  psychrometric  calculations  of  thermal  environment  begin,  usually, 
with  input  data  for  dry-  and  wet-bulb  temperature  of  moist  air.  These  two 
temperature  data  are  sufficient  to  describe  the  complete  thermodynamic 
state  of  the  moist  air  at  one  standard  atmosphere.  The  actual  calculations 
involving  heat-  and  mass-balance  within  a given  thermal  system  are,  however, 
more  readily  performed  with  the  dry-bulb  temperature  and  the  humidity  ratio, 
as  seen  in  equations  of  shelter  heat  balance.  The  method  used  during  this 
study  to  obtain  the  humidity  ratio  of  moist  air  from  its  dry-  and  wet -bulb 
data  is  described  as  follows: 

The  thermodynamic  wet-bulb  temperature  relation  2.7-7  can  be 
expanded  by  the  use  of  enthalpy  expression  2.7-6. 

—7—  [hs(t)  - ha(t)]  + ha(t)  + [Ws*(t')  - W]hw(t')  = hs*(t')  2.7-8 

vVt) 

In  the  above  expression,  Ws(t),  hs(t),  and  ha(t)  are  thermodynamic 
properties  at  dry-bulb  temperature  t,  while  Ws*(t/),  hs*(t/),  and  hw*(t/) 
are  thermodynamic  properties  evaluated  at  wet-bulb  temperature  t' . 

Rearranging  the  terms  in  equation  2.7-8,  the  humidity  ratio  W for 
dry-  and  wet-bulb  temperatures  t and  t/  can  be  expressed  as 

hs*(t')  - ha(t)  - hw*(t' )Ws*(t' ) 

W_  hs(t)  - ha(t)  - hw*(t')W8(t)  Ws(t)  2-7‘9 

The  calculation  of  the  thermodynamic  wet-bulb  temperature  from  a 
given  dry-bulb  temperature  and  humidity  ratio  is  also  possible  from 
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2.7-8  by  an  iterative  technique.  The  iterative  technique  found  successful 
during  the  course  of  this  investigation  is  the  inverse  interpolation 
formula  [6]  of  Newton  applied  to  Goff  and  Gratch  tables. 

2.8.  Effective  temperature  of  shelter  air. 

For  some  of  the  shelter  analyses,  effective  temperatures  have  been 
calculated  from  dry-  and  wet-bulb  temperatures  assuming  an  air  velocity  of 
less  than  20  fpm.  The  effective  temperature  chart  of  the  ASHVE*  has  been 
stored  in  the  computer  memory  in  a tabular  format,  and  a table -searching 
and  interpolative  subroutine  used  to  calculate  the  effective  temperature. 

3.  DESCRIPTION  OF  COMPUTER  PROGRAMS 

3.1.  Computer  models. 

Basically  four  different  computer  programs  have  been  developed  during 
this  study.  All  the  programs,  however,  essentially  employ  the  time-iteration 
technique  for  solving  transient  heat  conduction  equations  and  they  are 
designated  as  follows: 

M-(l)  Three-dimensional  model  with  composite  walls  and  with  separate 

initial  temperature  pattern  normal  to  the  six  boundary  surfaces. 

M-(2)  Three-dimensional  model  with  homogeneous  heat  conduction 

medium  having  separate  initial  temperature  patterns  normal  to 
the  six  boundary  surfaces . 

M-(3)  One -dimensional  compound  heat  conduction  model  for  composite 
regions  treating  the  six  exposures  separately. 

M-(4)  One -dimensional  compound  heat  conduction  model  for  six 
composite  regions  assuming  an  adiabatic  roof. 

Each  program  has  advantages  and  disadvantages,  as  discussed  in  the 
following  pages. 

* ASHVE  Guide  1950  Chapter  6 
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M-(l).  This  program  has  been  described  in  reference  [2],  and  was 
used  to  evaluate  the  NBS  family  shelter.  The  earth  temperature  field  sur- 
rounding the  shelter  was  divided  into  6 blocks,  such  as  shown  in  figure  1. 
This  block  system  enabled  the  program  to  account  for  situations  in  which 
some  wall  region  (s)  may  be  considerably  different  from  the  others  in  heat 
transfer  properties  and  earth  temperature. 

The  initial  temperature  in  each  block  was  programmed  only  in  the 
direction  normal  to  the  wall  surface,  however,  because  the  temperature 
profiles  were  usually  known  only  along  those  directions. 

Temperature  calculations  for  the  entire  earth  region  surrounding  the 
shelter  were  made  with  time-iterative  techniques  on  three-dimensional 
finite  difference  equations.  The  concrete  wall  (including  floor  and  ceiling) 
temperatures  were  calculated  by  one -dimensional  finite  difference  equations 
separately  for  each  wall,  assuming  that  lateral  temperature  variation  on 
the  interior  surface  for  a given  wall  could  be  neglected.  Each  inner 
wall-surface  temperature  was  determined  by  surface  heat  balance  equation 
2.2*  including  vapor  condensation  but  excluding  condensate  re-evaporation. 

The  shelter  psychrometric  condition,  dry-bulb  temperature,  dew-point  tem- 
perature, and  relative  humidity  were  then  evaluated,  based  upon  the  total 
heat  balance  equation  2.1. 

M-(2).  In  this  model,  thermal  properties  or  heat  transfer  charac- 
teristics around  the  shelter  air  space  all  were  assumed  homogeneous.  In 
other  words,  no  distinction  in  thermal  properties  was  made  from  the  concrete 
wall  to  the  soil,  or  from  one  wall  region  to  the  other,  as  in  Model  M-(l). 

The  earth  temperature  initialization  was  performed  only  in  the 
direction  normal  to  the  earth  surface.  The  finite  difference  scheme 
employed  for  the  three-dimensional  time-iteration  solution  of  the  heat 
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DIAGRAM  OF  HEAT  TRANSFER  MATRIX 
AROUND  FAMILY  SHELTER 
(SCHEMATIC) 


Fig*  1 


Schematic  diagram  of  the  matrix  used  for  computer  program  M-  (1)  . 
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conduction  equation  is  shown  in  figure  2.  The  homogeneity  assumption  of 
the  heat  conduction  equation  mentioned  makes  it  possible  to  analyze  only 
a one-quarter  segment  of  the  shelter,  because  of  the  symmetric  nature  of 
the  entire  system.  The  details  of  this  program  have  been  described  in 
reference  [1].  This  program  was  applied  to  the  NBS  6-man  family  shelter, 
the  Summerlin  shelter,  the  Broyles  shelter,  the  Napier  shelter,  and  the 
Reading  shelter,  for  the  purpose  of  comparing  the  calculated  shelter 
thermal  environment  with  the  experimentally  observed  data. 

M-(3)  . The  three-dimensional  effect  on  earth  heat  conduction 
around  the  shelter  becomes  less  and  less  significant  as  the  shelter  size 
increases.  A one -dimens ional  compound  system  was  developed  primarily 
for  large  shelters,  where  a major  portion  of  the  heat  flow  is  always 
normal  to  the  shelter  walls,  ceiling,  and  floor.  This  program  is  identi- 
cal with  M-(l),  except  that  the  corner  region  of  earth  and  concrete  is 
ignored,  and  the  one -dimens ional  finite  difference  equation  was  used  for 
earth  temperature  determinations . Comparison  with  the  two  previous 
programs  indicates  that  the  earth  temperature  computation  scheme  was 
drastically  simplified  in  the  model.  The  program  was  applied  to  the 
Summerlin  and  Reading  shelters. 

M-(4).  This  program  was  a modification  of  M-(3),  to  simulate 
basement  type  shelters  such  as  the  two  at  Ft.  Belvoir.  The 
heat  transfer  in  the  ceiling  region  of  basement  shelters  will  be  much 
less  significant  than  that  in  wall  and  floor  regions.  Thus,  in  this 
program,  the  outer  face  of  the  ceiling  layer  was  assumed  adiabatic. 
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Schematic  diagram  of  the  heat  conduction  region  used  for  computer 
program  M-  (2) . 


Fig.  2 
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The  calculation  of  the  psychrometr ic  condition  of  shelter  air  for 
this  program  was  designed  so  that  dry-  and  wet -bulb  and  dew-point  tempera- 
ture, as  well  as  shelter  effective  temperature,  can  be  computed  with  or 
without  the  air  conditioning  system  turned  on.  The  details  of  the  air 
conditioning  calculations  used  in  this  program  are  described  in  the 
appendix.  This  program  was  applied  to  200-man  and  1000-man  shelters  of 
Ft.  Belvoir,  which  were  tested  by  the  University  of  Florida. 

3.2.  Initialization  of  earth  temperature  surrounding  a protective 
shelter . 

Accurate  heat  transfer  calculations  for  the  early  part  of  an 
occupancy  period  for  these  shelters  are  extremely  difficult,  because  of 
the  uncertainty  about  the  earth  temperature  distribution  around  the 
shelter.  Particularly  for  the  shelters  that  had  been  installed  with  only 
a shallow  earth  cover,  the  temperature  variation  from  roof  region  to  the 
floor  region,  from  one  wall  region  to  another,  and  from  corner  to  flat 
surface  region  was  quite  appreciable,  and  very  complex  to  approximate 
mathematically . 

In  order  to  simulate  this  complex  and  three-dimensional  pattern  of 
the  initial  earth  temperature  around  the  shelter,  assuming  that  such 
three-dimensional  patterns  are  important,  the  earth  temperature  program 
becomes  highly,  and  perhaps  unnecessarily,  complicated.  Therefore,  during 
the  study,  a simplification  was  made  by  selecting  three  initial  tem- 
perature patterns,  as  follows: 

-1.  Earth  temperature  gradients  are  always  in  a direction 

normal  to  four  walls,  floor,  and  roof  surfaces  - block  mode. 
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1-2.  Predominant  earth  temperature  gradient  is  always  along  a 
vertical  path,  from  the  surface  downward  - vertical  mode. 

1-3.  Earth  temperature  is  virtually  constant  all  around  the 
shelter  - constant  mode. 

These  three  modes  of  the  initial  earth  temperature  patterns  have 
been  employed  during  the  calculation,  and  presented  in  this  report. 

However,  it  is  important  to  remember  that  the  detailed  earth  tem- 
perature profile  data  surrounding  shelters  seldom  will  be  available, 
even  when  it  is  important,  for  the  majority  of  actual  shelters.  For 
most  of  the  thermal  environment  calculations  of  underground  structures, 
it  is  usually  assumed  that  the  earth  has  a single  value  of  temperature. 
This  uncertainty  in  the  initial  earth  temperatures  is  one  of  the 
sources  of  error  in  predicting  the  thermal  environment  of  shelters, 
particularly  for  the  early  period  of  shelter  occupancy. 

4.  DESCRIPTION  OF  PROTOTYPE  SHELTERS 

Brief  descriptions  of  the  prototype  shelters  used  for  this  analysi 
and  their  characteristics  are  given.  If  the  mathematical  simulation  of 
shelter  thermal  environment  is  to  be  most  effective,  it  is  necessary  to 
secure  accurate  information  regarding  many  parameters  related  to  struc- 
tural, site,  climatic,  and  operational  characteristics  of  shelters. 

4.1.  NBS  shelter  . 

The  family  size  shelter  in  Washington,  D.  C.,  tested  by  the 
National  Bureau  of  Standards,  was  constructed  according  to  Bulletin  MP-15 
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of  the  Office  of  Civil  Defense  [3]  with  some  small  modifications.  It  was 
a concrete  wall  shelter  with  external  dimensions  of  12  feet  long,  9 feet 
4 inches  wide,  and  7 feet  6 inches  high  (not  including  hatch),  placed  in 
an  excavation  and  covered  with  2 feet  3 inches  of  earth.  The  interior 
dimensions  of  the  shelter  area  were  10  feet  8 inches  long,  8 feet  wide, 
and  6 feet  6 inches  high.  A 2-foot-wide  hatchway  was  installed  on  the 
north  side  of  the  occupancy  area  and  separated  from  the  8-  by  8-foot  living 
space  by  an  8-inch-thick  concrete  shielding  wall.  The  earth  surface  over 
the  shelter  was  grass  covered  and  partially  in  the  shade  of  neighboring 
trees.  The  soil  around  the  shelter  was  mostly  loam  and  clay;  its  density 
averaged  about  109  lb/ft3  , and  its  moisture  content  averaged  about  15 
percent  (dry  weight  basis). 

Ventilation  air  was  controlled  in  a neighboring  equipment  house  to 
specified  conditions  and  ducted  into  the  shelter  through  an  inlet  at  the 
mid -height  of  the  shielding  wall  facing  the  occupied  area,  and  the  shelter 
air  was  exhausted  through  an  outlet  located  on  the  opposite  wall  of  the 
hatchway . 

Six  simulated  occupants  (SIMOC)  were  carefully  designed  and  con- 
structed to  produce  sensible  and  latent  portions  of  metabolic  heat  as 
functions  of  shelter  area  temperature  (details  of  SIMOC  are  found  in 
reference  [3]).  Four  SIMOC 's  had  a nominal  heat  output  of  400  Btu/hr, 
and  the  other  two  had  heat  outputs  of  200  and  600  Btu/hr,  respectively. 

Complete  psychrometric  observations  were  made  at  the  ventilation 
air  inlet,  the  exhaust  air  outlet,  and  at  the  5-foot  level  above  the 
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geometrical  center  of  the  shelter  floor.  Temperature  observations 
obtained  during  this  test  included  earth  temperatures  around  the 
shelter  to  a maximum  distance  of  4 feet  away  from  the  exterior  surface  of 
the  concrete  wall  for  all  four  walls  and  the  floor.  In  the  shelter  roof 
region,  earth  temperature  was  studied  with  four  thermocouples  located 
6 inches,  12  inches,  18  inches,  and  24  inches  from  the  external  surface 
of  the  shelter  ceiling. 

Five  tests  were  conducted  with  variations  in  test  duration,  venti- 
lation rate,  number  of  occupants,  ventilation  air  condition,  and  sur- 
rounding earth  temperature.  Test  1 was  conducted  with  no  occupants,  with 
42  cfm  of  ventilation  air.  Tests  2,  3,  and  4 were  each  conducted  with 
six  SIMOC's,  and  with  ventilation  air  rates  of  0,  18,  and  42  cfm,  respec- 
tively. Test  5 was  undertaken  to  simulate  winter  conditions  of  occupancy, 
ventilation  air,  and  earth  temperature,  using  six  SIMOC's  and  18  cfm  venti- 
lation air  . 

4.2.  Summerlin  shelter  [4]. 

This  shelter  was  a welded  steel  structure  located  entirely  below  the 
finished  grade  line,  with  a 30-inch  earth  cover  over  the  roof,  in  a rural 
area  of  Gainesville,  Fla.,  and  thermally  isolated  from  other  buildings. 

The  occupancy  area  dimensions  were  28  feet  7 inches  long,  and  7 feet  8 
inches  wide.  The  shelter  roof  was  arched  over  the  wall,  and  the  maximum 
ceiling  height  at  the  middle  of  the  arch  was  7 feet  1\  inches.  A 4-  by 
3-foot  hatchway  was  located  at  one  end  of  the  shelter.  The  earth  around 
the  shelter  was  a mixture  of  sand  and  loam,  and  the  earth  cover  was  bare 
at  the  time  of  the  shelter  environmental  tests.  The  moisture  content  (dry 
weight  basis)  of  the  soil  samples  analyzed  ranged  from  9 percent  to  19 
percent;  their  dry  density  ranged  from  103  lb/ft3  to  111  lb/ft°  . 
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Ventilation  air  conditions  were  generated  in  an  equipment  trailer 
located  outside  the  shelter,  and  supplied  to  the  shelter  by  a flexible 
tube  passing  through  the  entrance  hatch.  The  ventilation  air  inlet  was 
near  the  north  end  of  the  shelter,  and  shelter  air  was  exhausted  from  a 
stack  located  at  the  southeast  corner. 

Two  series  of  environmental  tests  were  conducted,  using  18  SIMOC's 
(ref.  3)  of  NBS  400  Btu/hr  type;  one  during  July  1962  with  ventilation 
air  of  200  cfm  and  typical  Florida  summer  outdoor  psychrometr ic  condi- 
tions; and  the  other  during  April  1963  with  ventilation  rates  of  54  and 
216  cfm  at  August  and  April  psychrometr ic  Conditions.  The  simulated  August 
psychrometric  conditions  used  for  the  tests  were  diurnal  cycles  of  96  °F 
maximum  dry-bulb  temperature,  80  °F  minimum  dry -bulb  temperature,  and 
79.8  °F  average  dewpoint  temperature.  The  simulated  April  psychrometric 
conditions  used  for  the  test  were  diurnal  cycles  of  76  °F  maximum  dry- 
bulb  temperature,  61  °F  minimum  dry -bulb  temperature,  and  59.5  °F  average 

dewpoint  temperature . 

4.3.  The  Broyles  shelter  [ 4]  8 

The  Broyles  Shelter,  at  Gainesville,  Fla.,  was  so  constructed  that 
a portion  of  the  shelter  was  below  grade,  and  3 feet  of  earth  was  mounded 
around  the  above -grade  portion.  The  interior  of  the  shelter  had  a ceiling 
height  of  7.33  feet,  and  the  floor  area  measured  16  by  7.75  feet,  for  a 
total  area  of  124  square  feet.  The  floor,  placed  on  a plastic  membrane, 
was  of  waterproof  concrete  reinforced  with  1/2-inch  steel  rods  on  12-inch 
centers.  The  walls  were  hollow  concrete  blocks  with  conventional  mortar 
joints,  and  the  cavities  were  filled  with  waterproof  cement  as  the  walls 
were  constructed.  The  roof  was  a concrete  slab  reinforced  with  3/8-inch 
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steel  rods  on  12-inch  centers.  The  concrete  was  treated  with  a water- 
proofing material  at  the  time  it  was  mixed.  This  shelter  was  shaded  by 
surrounding  trees,  and  a heavy  layer  of  sod  and  green  grass  covered  it  and 
the  surrounding  ground.  The  surrounding  earth  was  a mixture  of  sand  and 
loam,  whose  dry  density  ranged  from  95  lb/ftJ  to  104  lb/ft3  , and  moisture 
content  from  1 percent  to  12.5  percent. 

For  most  of  the  thermal  environment  test  period,  12  SIMOC's  (ref.  3) 
of  NBS  400  Btu/hr  type  were  employed.  The  ventilation  air  was  processed 
outside  the  shelter  to  simulate  diurnal  cyclic  conditions  of  96  °F  maxi- 
mum dry-bulb  temperature,  80  °F  minimum  dry-bulb  temperature,  and  70  °F 
average  dew-point  temperature,  and  was  forced  into  the  shelter  by  a fan 
through  an  air  supply  duct  that  was  laid  along  the  inside  surface  of  the 
longitudinal  shelter  wall,  and  which  had  four  equally  spaced  outlets.  The 
shelter  air  was  exhausted  through  two  outlets  located  at  opposite  ends  of 
the  shelter. 

In  addition  to  the  regular  psychrometr ic  measurements  of  inlet 
ventilation  air,  shelter  exhaust  air,  and  shelter  air  at  the  geometrical 
center,  several  measurements  were  made  of  the  shelter  inner  wall  surface 
temperature  and  surrounding  earth  temperature . 

The  shelter  was  tested  in  four  successive  phases: 

Phase  1:  A ventilation  rate  of  3 cfm  per  person  was  supplied, 

and  a fan-and-coil  unit  simultaneously  cooled  the  recirculated  air  by 
4 gpm  of  well  water  at  an  inlet  water  temperature  of  71.5  °F.  During 

the  first  two  days,  sensible  and  latent  heat  was  removed  from  the  shelter  by 
this  fan-and-coil  unit  at  a rate  of  6440  Btu/hr*  This  cooling  capacity 
exceeded  the  total  heat  supplied  to  the  SIMOC^s  by  4800  Btu/hr* 
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Phase  2.  The  well-water  coil  was  cut  off  and  10  cfm  per 

person  of  ventilation  air  was  supplied  for  a period  of  12  days. 

Phases  3 and  4.  These  two  phases  were  operated  under  the 
same  conditions  as  phase  2,  except  that  the  ventilation  rate  in  phase  3 
was  6 cfm  per  person  (for  48  hours)  and  in  phase  4 was  3 cfm  per  person 
(for  48  hours ) . 

4.4.  Napier  shelter  [4]  . 

The  Napier  shelter  was  a 100-occupant  community  shelter  designed 
for  a group  of  residents  in  a subdivision  adjoining  the  city  of 
Gainesville,  Fla.  The  floor  level  was  5 feet  below  grade  and  30  inches 
of  earth  covered  the  roof.  The  floor  slab  was  wire-mesh  reinforced, 
4-inch-thick  concrete  poured  over  a waterproof  plastic  membrane.  The 
walls  were  of  8-inch-thick  hollow  concrete  blocks,  whose  cavities  had 
reinforcing  rods  placed  vertically  through  then  at  selected  intervals 
and  were  then  filled  with  concrete.  Outside  dimensions  were  20  feet  wide 
and  85  feet  long.  The  interior  floor  area  was  1561  square  ffeet . Rein- 
forced prestressed  concrete  T-beams  placed  on  top  of  the  shelter  walls, 
each  in  contact  with  the  beams  parallel  to  it,  formed  its  roof.  The 
surrounding  earth  was  mostly  clay,  whose  dry  density  varied  from  76  lb/ft3 
to  109  lb/ft3  , while  the  moisture  content  ranged  from  3.4  percent  to 
34.2  percent  (dry  weight  basis).  The  earth  surface  over  the  shelter  was 
bare  with  several  patches  of  weeds. 

For  this  test,  100  SIMOC’s  (ref.  3)  of  NBS  400  Btu/hr  type  were 
employed.  The  ventilation  air  was  conditioned  in  the  equipment  trailer 
outside  the  shelter  to  represent  a diurnal  cycle  of  a Florida  summer: 

96  °F  maximum  dry-bulb  temperature,  80  ®F  minimum  dry-bulb  temperature, 
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and  73  °F  average  dewpoint  temperature.  It  was  forced  into  the  shelter 
at  two  stations  at  the  ceiling  level  in  the  south  wall,  approximately  8 
feet  from  the  southwest  and  southeast  corners  of  the  shelter.  The 
ventilation  air  was  discharged  in  the  direction  of  the  north  wall;  and 
after  passing  through  the  main  chamber  was  exhausted  through  the  entry 
doorway  and  entry  chamber.  Regular  psychrometric  measurements  were  made 
of  shelter  air,  ventilation  air,  and  exhaust  air.  In  order  to  study  the 
three-dimensional  earth  temperature  profile  during  the  simulated  occu- 
pancy period  of  the  shelter,  numerous  measurements  of  earth  temperature 
were  made  around  the  northwest  and  southwest  corner  regions. 

The  test  was  conducted  in  six  phases : 

Phase  1.  Ventilation  rate  was  maintained  at  3 cfm  per  person 
for  48  hours,  during  which  a well-water  coil  using  12  gpm  of  72.2  °F  water 
was  in  operation.  The  measured  total  cooling  capacity  of  the  well-water 
coil  during  the  phase  averaged  20,700  Btu/hr,  amounting  to  51.6  percent 
of  the  total  heat  released  by  SIMOC's. 

Phase  2 . Ventilation  air  rate  was  increased  to  6 cfm  per 
person  and  the  well  water  coil  was  shut  off  for  this  period  of  97  hours. 

Phase  3.  With  all  other  conditions  being  identical  with  those 
of  phase  2,  30  NBS  SIMOC's  were  replaced  by  one  MASS  SIMOC  of  MRD  [4]  for 
42  hours.  (One  MASS  SIMOC  has  an  adjustable  output  from  1-40  SIMOC's  of 

NBS  400  Btu/hr  type.) 

Phase  4.  The  MASS  SIMOC  was  replaced  by  30  NBS  SIMOC's  for 

51  hours . 

Phase  5.  Ventilation  air  rate  was  increased  to  8.05  cfm  per 


person  with  100  NBS  SIMOC's  for  116  hours. 
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Phase  6,  The  total  ventilation  rate  remained  at  805  cfm 
and  100  NBS  SIMOC’s  were  used  together  with  1 MASS  SIMOC  for  36  hours, 
representing  a total  of  140  occupants. 

4.5.  Reading  shelter  [4]. 

This  was  a community  shelter  located  in  a park  owned  by  the  city 
of  Reading,  Pa.  It  was  constructed  in  a hillside  to  take  advantage  of  a 
thick  earth  cover,  and  was  basically  a rectangular  parallepiped  56  feet 
long,  17  feet  4 inches  wide,  with  a ceiling  height  of  7 feet  8 inches. 

The  ends  of  the  shelter  were  connected  to  separate  tunnel-like  entry 
corridors  so  that  two  right  angle  turns  were  formed  in  each  of  these 
passageways.  The  floor  was  4 feet  6 inches  below  the  grade  level  that 
existed  prior  to  construction,  and  the  roof  and  all  sides  had  a minimum 
earth  covering  of  JO  inches.  All  bearing  walls,  the  roof,  and  the 
floor  were  constructed  of  concrete  reinforced  with  steel  bars.  Water- 
proofing was  applied  to  the  external  surfaces  of  the  shelter  during 
construction.  There  were  many  internal  partition  walls  constructed  of 
hollow  concrete  blocks  with  mortar  filled  voids,  for  rooms  of  various 
purposes,  such  as  storage,  first  aid,  mechanical  equipment,  and  lavatory. 
The  surrounding  earth  was  sandy  loam  of  dry  density  between  89  lb/ft3  and 
121  lb/ft°,  and  moisture  content  between  10  percent  and  22  percent  (dry 
weight  basis).  The  earth  surface  at  the  time  of  testing  was  grass  under 
snow.  Ventilation  air  artificially  created  in  an  equipment  trailer 
outside  the  shelter  was  carried  to  the  shelter  through  a duct  connected 
to  a stack  that  under  normal  shelter  operation  would  be  utilized  as  an 
exhaust  stack.  A temporary  distribution  duct  for  ventilation  air  supply 
was  installed  at  the  ceiling  level  along  the  west  wall  of  the  shelter. 
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Three  air  outlets  were  installed  in  the  temporary  duct  at  equally  spaced 
intervals.  The  air  was  exhausted  through  a stack  located  in  the  center 
of  the  ceiling  of  the  equipment  room,  which  was  almost  longitudinally  at 
the  opposite  end  from  the  air  supply  duct. 

The  ventilation  air  conditions  used  for  this  test  varied  considerably; 
however,  the  dry-bulb  temperature  was  generally  maintained  between  32  and 
40  °F,  with  an  average  dewpoint  temperature  of  JO  °F.  The  simulated 
occupants  used  were  two  MASS  SIMOC’s  (ref  4),  which  were  capable  of 
producing  sensible  and  latent  metabolic  heat  equivalent  to  120  sedentary 
adults . 

Psychrometric  observation  stations  were  located  at  the  ventilation 
air  inlet,  geometric  center  of  the  shelter,  and  exhaust  air  outlet. 
Thermocouples  were  used  to  measure  interior  surface  temperatures  and 
several  ground  temperatures  extending  downward  45  inches  distance  from  the 
top  of  the  floor  slab  and  horizontally  from  the  outside  surface  of  the  west 
wall.  The  roof  region  ground  temperature  was  also  measured  at  several 
distances  from  the  outside  surface. 

The  test  was  divided  into  seven  phases,  described  as  follows: 

Phase  1.  Ventilation  at  the  rate  of  150  cfm  was  supplied 
to  the  shelter  in  the  manner  described  previously,  and  the  two  MASS  SIMOC's 
were  adjusted  to  deliver  heat  and  moisture  equivalent  to  a total  of  50 
sedentary  adults.  This  phase  continued  from  February  26  to  March  4,  1963. 

Phase  2.  Ventilation  rate  was  then  decreased  to  75  cfm,  or 
1.5  cfm  per  person,  with  operating  conditions  identical  to  those  of 
phase  1.  This  phase  lasted  approximately  5 days. 


- 24  - 


Phase  3.  The  mass  SIMOC’s  were  readjusted  to  produce  total 
equivalent  metabolic  heat  for  100  occupants.  The  total  ventilation  air 
rate  was  150  cfm  during  this  phase  of  study,  lasting  three  days. 

Phase  4.  During  the  three-day  period  that  followed  phase  3, 
from  March  12  to  14,  shelter  occupancy  was  50  simulated  occupants.  The 

ventilation  air  at  a rate  of  83I  cfm  was  supplied  by  the  blower  that 
originally  had  been  installed  in  the  shelter  as  a part  of  the  permanent 

facilities.  The  shelter* s own  air  handling  systems  provided  the  distri- 
bution and  exhaust  of  the  shelter  air.  The  purpose  of  this  phase  of 
the  test  was  to  determine  if  the  original  equipment  for  supplying  air, 
and  the  air  distribution  and  venting  systems  were  adequate. 

Phase  5 . During  this  phase,  ventilation  air  at  the  rate  of 
718  cfm  was  supplied  for  50  simulated  occupants . The  air  distribution 
and  exhaust  systems  were  changed  during  this  phase,  details  of  which  are 
not  important  for  the  purpose  of  this  report. 

Phases  6 & 7 . These  phases  were  conducted  for  sealed-up 
conditions  without  ventilation  for  50  simulated  occupants  for  24  hours, 
and  100  cimulated  occupants  for  58  hours,  respectively. 

4.6.  Ft.  Belyoir  200-man  shelter. 

This  was  an  experimental  200-man  shelter  designed  and  built  by  the 
Protective  Structures  Development  Center  at  Ft.  Belvoir,  Va . it  was  a 
two-story  reinforced  concrete  structure  with  one  story  below  ground  level. 
Because  the  building  had  been  occupied  by  personnel  of  the  Protective 
Structures  Development  Center  prior  to  the  environmental  testing,  the 
temperature  and  humidity  had  been  comfort  conditioned.  Only  the  base- 
ment area  was  used  for  the  shelter  environment  test,  while  the  upper 
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floor  area  was  heated  to  simulate  an  adiabatic  roof  condition.  The 
shelter  walls  were  10  inches  thick,  the  floor  slab  6 inches  thick,  and 
the  ceiling-floor  slab  8 inches  thick.  There  was  a 12 -inch-thick  fill 
of  gravel  between  the  external  surface  of  concrete  of  the  basement  and 
the  earth  on  all  sides,  except  the  north  end.  The  earth  around  the 
shelter  was  a mixture  of  sand  and  clay,  its  dry  density  varying  from 
85  lb/ft3  to  120  lb/ft3  and  its  moisture  content  from  9 percent  to  23 
percent  (dry  weight  basis).  The  earth  surface  was  composed  of  sod  on 
all  sides  of  the  shelter,  except  for  a portion  of  the  north  end  where 
a bituminous  concrete  driveway  was  located.  Interior  dimensions  were 
37  feet  2 inches  by  37  feet  2 inches,  with  a ceiling  height  of  9 feet  6 
inches.  Excluding  a first  aid  room  and  stairwell,  a net  usable  floor 
area  of  1032  square  feet  was  available  for  100  simulated  occupants. 

There  were  two  family-type  basement  shelters  that  were  built  for  display 
purposes  in  the  test  room,  one  of  sand-filled  concrete  block  walls  and 
the  other  a triangular-shaped  wooden  lean-to  filled  with  sand  and  sand 
bags.  They  may  have  had  a considerable  effect  in  absorbing  heat  during 
the  first  few  days  of  the  environmental  test. 

The  temperature-  and  humidity-controlled  ventilation  air  (dry-bulb 
cycles  93  °F  ^ 76  °F  and  wet-bulb  78  °F~  74  °F  were  in  phase  with 
maximum  at  2 p.m.  and  minimum  at  4 a.m.)  was  produced  in  the  equipment 
trailer  parked  outside  the  shelter  and  introduced  into  the  shelter  from 
ceiling  level  near  the  center  of  the  occupied  space.  The  exhaust  outlet 
was  also  in  the  ceiling,  approximately  12  feet  from  the  air  inlet.  The 
thermal  environment  of  the  lower  shelter  space  was  measured  with  100 
SIMOC^s  of  the  NBS  type  (ref.  3)* 
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The  test  was  conducted  in  five  phases,  in  which  the  ventilation 
air  rates  were  3 cfm  per  person  for  the  first  phase  (6  days),  6 cfm  per 
person  for  the  second  phase  (2  days),  9 cfm  per  person  for  the  third 
phase  (3  days),  18  cfm  per  person  for  the  fourth  phase  (2  days),  and  27 
cfm  per  person  for  the  fifth  phase  (4^  days). 

Psychrometric  measurements  were  made  at  the  ventilation  air 
inlet,  air  exhaust,  east  geometric  center,  and  west  geometric  center  of 
the  shelter.  Temperature  measurements  were  obtained  of  various  inner 
surfaces,  together  with  surrounding  earth  temperatures  to  a distance  of 
4 feet  from  the  external  surface  of  the  concrete  structure . 

4.7.  Ft.  Belvoir  1000-man  shelter. 

This  shelter  was  built  by  the  Protective  Structures  Development 
Center  at  Ft.  Belvoir,  Va . It  was  a two-story  reinforced  concrete 
building  with  one  story  below  ground  level.  Prior  to  testing  the  building 
had  been  occupied  by  personnel  of  the  Protective  Structures  Development 
Center,  and  its  thermal  environment  had  been  controlled  at  a comfort  air 
condition  by  a central  heating  and  cooling  plant.  Only  the  basement 
area  was  used  for  test,  while  the  upper  floor  level  was  heated  to  simu- 
late an  adiabatic  roof  condition.  The  usable  floor  area  of  the  basement 
was  5400  and  was  given  a simulated  occupancy  of  11  MASS  SIMOC’s  (ref. 
4).  The  shelter  walls  were  10  inches  thick,  the  floor  6 inches  thick,  and 
the  ceiling-f loor  slab  8\  inches  thick.  The  interior  dimensions  were  74 
feet  4 inches  by  74  feet  4 inches,  with  a ceiling  height  of  9 feet  6 
inches.  The  surrounding  earth  was  sandy  clay  and  ranged  in  moisture 
content  from  10  to  20  percent  (dry  weight  basis),  and  in  dry  density 
from  86  lb/ft3  to  104  lb/ft3  depending  upon  the  depth,  as  well  as  upon 
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the  location  around  the  shelter.  During  the  backfilling  and  mounding 
operations,  a 12-inch-thick  fill  of  gravel  was  placed  adjacent  to  all 
the  basement  walls. 

The  ventilation  air  was  conditioned  by  the  equipment 
trailer  to  yield  a local  design  diurnal  cycle  of  dry-bulb  temperature, 

94  °F  maximum,  76  °F  minimum,  and  wet-bulb  temperature,  78  °F  maximum, 

74  °F  minimum.  The  maximum  and  minimum  conditions  of  dry-  and  wet -bulb 
temperatures  were  both  at  2 p.m.  and  4 a.m.,  respectively.  The  ventila- 
tion air  inlet  and  exhaust  air  outlet  were  both  near  a corner  of  the 
space  and  separated  from  each  other  by  approximately  12  feet  for  most  of 
the  test  period. 

The  test  was  divided  into  four  major  phases,  excluding  an  initial 
short  period  of  purging.  Phase  1 was  a sealed-up-condition  test  of 
approximately  8 hours  without  ventilation,  whereas  the  ventilation  rates 
for  remaining  phases  were  varied  from  16.8  efm  per  person  for  34  days, 
to  5 efm  per  person  for  1 day.  The  shelter  thermal  environment  was 
measured  by  dry-  and  wet-bulb  thermometers  at  the  exhaust  duct,  the 
northeast  area,  the  southwest  area,  and  the  geometric  center  of  the 
shelter,  and  numerous  thermocouple  temperature  readings  were  taken  for 
many  parts  of  the  wall  surface  and  the  earth  region  extending  outward 
normally  to  the  walls  and  floor  to  a distance  of  4 feet  from  the 
exterior  surfaces. 

5.  SHELTER  CALCULATIONS 

As  discussed  in  section  2 of  this  report,  the  calculated  shelter 
thermal  environment  is  affected  by  several  parameters,  such  as  thermal 
conductivity  and  thermal  diffusivity  of  surrounding  walls  and  earth. 


- 28 


and  heat  transfer  coefficients  along  the  inner  surfaces  of  the  shelter  and 
the  ground  surface.  Based  upon  the  characteristics  of  prptotype  shelters 
and  their  test  conditions  described  in  the  previous  section,  table  1 is  pre- 
pared to  summarize  the  parameters  used  for  the  shelter  calculations.  However 
many  of  the  heat  transfer  parameters  selected  for  the  computation  cannot  be 
too  precise.  The  following  considerations  were  given  for  assigning  numerical 
values  to  the  parameters  listed  in  table  1. 

5.1.  Physical  dimensions  (size  and  heat  transfer  area). 

The  size  of  the  shelter  is  expressed  in  the  overall  internal 
dimensions,  which  are  the  overall  external  dimensions  less  the  thickness 
of  walls,  ceiling,  and  floor.  The  heat  transfer  area  used  for  the  com- 
putation was  calculated  from  the  internal  dimension  only,  thus  ignoring 
the  complex  pattern  of  partitioning  and,  consequently,  the  heat  absorp- 
tion by  partition  walls  or  columns. 

5.2.  Simulated  occupants. 

The  simulated  occupants  (SIMOC's)  used  for  the  experiment  were 
designed  to  simulate  the  heat  output  of  human  bodies  and  each  generated  a 
total  heat  of  400  Btu/hr,  regardless  of  temperature.  The  sensible  heat 
component  was  regulated  according  to  the  shelter  air  temperature  in  the 
manner  described  in  section  2.  The  number  of  active  SIMOC '.s  was  varied 
during  some  tests  (Napier  and  Reading  shelters)  by  turning  the  power 
off  and  on  to  a part  of  them.  For  this  reason,  the  area  per  person 
data  shown  in  table  1 for  these  tests  are  not  constant. 

In  order  to  simulate  the  change  in  number  of  SIMOC's  during  some 
of  the  tests,  the  NBS  computer  program  was  constructed  to  take  the  number 
of  occupants  as  a time  variable  instead  of  a constant. 


Table  1.  UNDERGROUND  FALLOUT  SHELTER  CALCULATIONS 
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5.3.  Ventilation  air. 

Observed  hourly  data  on  dry-bulb  and  dewpoint  temperatures  or 
dry-bulb  and  wet-bulb  temperatures  of  the  supply  air  duct  were  used  as 
the  input  data,  either  in  the  2-hour  step  or  the  4-hour  step  calculations. 
Since  the  University  of  Florida  varied  the  ventilation  air  rate  during 
their  tests,  the  computer  program  was  designed  to  accept  the  ventilation 
rate  as  time  dependent  input  data  . These  variations  of  the  number  of 
SIMOC's  and  ventilation  rates  during  the  tests  are  the  reason  that  the 
per  capita  ventilation  rate  indicated  in  table  1 was  not  a constant  for 
some  of  the  tests . 

The  ventilation  air  temperature  and  humidity  conditions  were 
varied  with  respect  to  time,  closely  following  a prescribed  diurnal 
cyclic  pattern  in  most  cases.  The  entries  in  table  1 show  approximate 
ranges  of  the  temperature  levels  employed  during  the  test,  whereas  actual 

hourly  values  were  used  in  the  calculations. 

5.4,  Thermal  properties  of  earth  and  wall. 

For  all  of  the  prototype  shelters.  Soil  samples  were  taken  from 
several  representative  spots  around  the  shelter  at  selected  depths  and 
were  analyzed  with  respect  to  soil  classification,  dry  density,  and  water 
content  (dry  weight  basis);  their  typical  characteristics  are  shown  in 
table  1.  The  thermal  conductivities  and  dif fusivities  of  various  types 
of  soils  are  usually  presented  as  functions  of  moisture  content  [5,6]. 

Such  charts  were  consulted  in  arriving  at  the  values  of  thermal  conduc- 
tivity, k,  and  thermal  diffusivity,  or,  in  table  1. 

For  the  NBS  shelter,  the  thermal  diffusivity  value  of  0.026  ft2  /hr. 
for  tests  1 and  2 was  estimated  from  a phase  angle  shift  of  the  earth 
temperature  cycles  at  two  different  depths.  The  diffusivity  values  of 
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0.022  ft2 /hr  for  tests  3 and  4,  and  0.03  ft2 /hr  for  test  5 were  evaluated 
by  a numerical  technique  similar  to  that  reported  by  Beck  [7].  The  method 
is  basically  a reversed  use  of  finite  difference  solutions  of  the  heat 
conduction  equation  for  a semi-infinite  solid.  For  the  case  of  composite 
wall  models,  the  thermal  diffusivity  and  conductivity  of  concrete  were 
assumed  as  0.036  ft2 /hr  and  1.1  Btu/hr  ft2 (deg  F/ft),  respectively. 

5.5.  Surface  Heat  Transfer  Coefficients 

The  surface  heat  transfer  coefficient  consists  of  a radiative  por- 
tion and  a convective  portion.  For  the  heat  transfer  at  the  interior  shel- 
ter surfaces,  the  radiation  portion  plays  a predominant  role,  since  the  air 
velocity  over  the  shelter  inner  surface  usually  is  very  small.  In  addi- 
tion, the  surface  heat  transfer  coefficient  along  any  one  of  the  vertical 
walls  at  a given  time  would  vary  considerably  from  the  bottom  to  the  top, 
due  to  the  varying  nature  of  layer  pattern,  as  well  as  that  of  the  radia- 
tion heat  exchange  geometry  and  air  velocity,  and  also  because  of  the 
difference  in  local  temperature  distribution.  However,  it  is  probable 
that  the  local  variation  of  the  surface  heat  transfer  coefficient  along  a 
given  surface  may  be  of  the  same  order  of  magnitude  as  the  convective  heat 
transfer  coefficient  itself.  Several  attempts  were  made  during  this  study  to 
obtain  the  surface  heat  transfer  coefficients  for  the  experimental  obser- 
vations of  the  shelter  wall  heat  conduction.  Since  the  NBS  test  shelter 
was  equipped  with  a heat  flow  meter  at  the  geometric  center  of  each  of 
all  the  inner  surfaces,  and  since  the  inner  surface  temperatures  and 
shelter  air  temperatures  were  simultaneously  measured,  it  was  possible 
to  calculate  the  heat  transfer  coefficients.  However,  the  accuracy  of 
this  procedure  is  questionable,  because  (a)  the  heat  flow  meter  reading 
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at  the  geometric  center  of  an  inner  surface  did  not  necessarily  yield 
an  average  heat  transfer  coefficient  for  the  entire  surface,  because  of 
the  local  variation  of  heat  flow,  and  (b)  where  condensation  of  shelter 
air  moisture  was  taking  place  simultaneously,  the  heat  flow  meter  would 
read  the  total  heat  flux,,  which  is  not  proportional  to  the  temperature 
difference  between  the  surface  and  air.  In  fact,  these  difficulties 
were  realized  in  the  NBS  shelter,  since  the  ratio  agreement  between  the 
shelter  total  heat  conduction  estimated  by  the  heat  balance,  and  the 
heat  conduction  based  upon  the  heat  flow  meter  ranged  from  52  percent 
to  108  percent  in  various  tests  [2]. 

Nevertheless,  the  values  of  heat  transfer  coefficients  listed  for 
the  NBS  shelter  tests  1,  2,  3>  and  4,  and  the  Summerlin  shelter  test  we^e 
estimated  from  the  test  3 data  of  the  heat  flow  meter  readings,  adjusted 
by  the  vapor  transfer  due  to  the  difference  of  air  humidity  ratio  between 
the  air  and  the  wet  surface,  and  the  Lewis  relation  of  heat  and  mass  transfer. 

A low  value  for  the  floor  coefficient  was  observed.  This  result 
was  to  be  expected  since  the  floor  surface  was  typically  colder  than  the 
air  immediately  above  it  and  the  downward  connection  heat  transfer  rate  would 
be  very  low  under  these  conditions.  Moreover  the  effective  heat  transfer 
area  of  the  floor  was  considerably  reduced  by  the  presence  of  simulated, 
occupants.  The  same  values  for  the  surface  coefficient  were  applied  to 
the  test  5 condition  of  the  NBS  shelter,  but  the  agreement  between  the 
observed  and  calculated  inner  surface  temperatures  was  rather  poor 
when  this  low  value  of  the  surface  coefficient  was  used. 
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Several  other  combinations  of  the  surface  heat  transfer  coefficients 
for  NBS  test  5 were  tried  and  the  results  are  summarized  in  table  2,  which 
will  be  elaborated  later  in  this  section.  For  the  many  other  prototype 
shelters,  heat  flow  meter  readings  were  not  included,  or  were  not  usable, 
so  that  the  detailed  heat  transfer  analysis  on  the  surface  by  surface 
basis  was  discarded.  However,  the  average  heat  transfer  coefficients 
for  the  entire  inner  surface  were  obtained  from  the  sensible  heat 
balance  calculation  for  the  Summerlin  shelter  during  the  moderate  weather 
condition  test,  the  NBS  test  3>  the  NBS  test  5,  the  Ft.  Belvior  1000-man 
shelter,  and  the  Napier  shelter. 

As  indicated  in  figures  3 and  4,  the  result  of  this  analysis  on 
the  overall  sensible  inner  surface  heat  transfer  coefficient  shows  a 
greatly  fluctuating  pattern.  The  accuracy  of  the  calculated  values  of  the 
surface  heat  transfer  coefficient  is  inherently  related  to  the  accuracy  of 
measuring  air  and  surface  temperatures.  Since  radiation  would  usually  be 
present  and  wetted  surfaces  are  sometimes  involved^  the  error  in  temperature 
measurement  could  easily  be  a signigicant  part  of  the  observed  temperature 
difference  for  differences  of  2 °F  or  less.  A significant  trend  can  be 
observed,  however,  from  Figure  3 and  4 that  the  overall  inner  surface 
sensible  heat  transfer  coefficient  seems  to  increase  as  the  temperature 
difference  between  the  air  and  inner  surface  decreases.  The  majority  of  the 
sensible  heat  transfer  coefficients  are  in  the  neighborhood  of  1.0,  except 
for  the  Napier  shelter.  Nevertheless,  the  inner  surface  heat  transfer 
coefficients  selected  for  the  Broyles  and  Napier  shelters  of  table  1 are 
somewhat  arbitrary,  whereas  those  given  to  the  Reading  and  Ft.  Belvior 
shelters  are  estimated  from  ref,  8,  which  was  not  available  prior  to 
the  last  three  shelter  calculations.  The  possible  effects  of  the  various 


Table  2.  EFFECT  OF  INNER  SURFACE  HEAT  CONDUCTANCE 
FOR  NBS  SHELTER  TEST  5 


- 34  - 


o 

is 


o 

OOOOOO^-iOO 


o 

OOOOOOt-iOO 


u 


CM 


• 

i — 1 

CN 

00 

0 

P- 

CM 

CM 

in 

CO 

1 

O 

Mt 

00 

0 

in 

1 — 1 

CM 

co 

4-1 

0 

M3 

r-4 

CM 

00 

CM 

CO 

CM 

Mt 

Mt 

i ^ 

1 — l 

CO 

r*- 

r-4 

0 

T“ — 1 

M3 

Mt 

1+4 

T— 4 

cd 

p". 

n 

CM 

in 

CM 

0 

Mt 

CM 

in 

! O 

in 

CM 

r — 1 

CM 

r-» 

CO 

CM 

1 — 1 

CO 

0 

00 

r^. 

n- 

P~-- 

P- 

00 

i n* 

r>» 

p- 

p-- 

C'- 

f'- 

r-. 

O 

i 


© 

S3 

pi 


. o ' in 

CO  • V 

M3  0 y 1 00 

OO  ! 00 


I 


• 

p^ 

1 1 

00 

CM 

00 

M3 

M3 

co 

in 

! CM 

CO 

UO 

r-4 

00 

00 

co 

CO 

ON 

0 

0 

O 

CO 

CM 

CO 

CO 

I-"- 

CM 

UO 

i CM 

t-4 

CO 

00 

Mt 

UO 

r-- 

r-4 

O 

> 

T 1 

%r4 

O 

i 

cd 

vD 

00 

O 

T— 1 

0 

M3 

00 

O 

r-4 

00 

ON 

r-4 

r-4 

r-4 

O 

i-l 

CM 

03 

O 

in 

uo 

M3 

M3 

M3 

in 

UO 

M3 

M3 

1 UO 

UO 

M3 

M3 

M3 

UO 

M3 

M) 

M3 

P 

CM 

1 

j 

M-( 

4-J 

4-1 

4-4 

! u 

1 

•H 

1 

• | 

00 

co 

X) 

u 

1 id 

CO 

• 

• 

M3 

i 4-1 

M3  1 

CO 

— 

— 

— 

-> 

1 ^ 

~ 

T“H 

1 

O ' 

m 

M) 

id 

P 

! 

s 

•U 

1 

} 

5-1 

pq 

! 

a> 

• 

r-4 

M3 

r- 

M3 

CO 

00 

00 

Ml" 

Mt 

' UO 

M3 

M3 

00 

ON 

t-4 

00 

CO 

x: 

0 

a 

CM 

O 

00 

Mt 

Mt 

r-* 

r-'. 

O 

ON 

UO 

Mt 

00 

ON 

CM 

ON 

M3 

CO 

-u 

• 

r4 

uo 

cd 

in 

M3 

M3 

M3 

p-- 

<fr 

M3 

M3 

uo 

M3 

P-^ 

00 

r^- 

00 

M3 

00 

C- 

P^ 

O 

n 

UO 

UO 

UO 

uo 

uo 

UO 

uo 

UO 

UO 

uo 

UO 

in 

UO 

in 

UO 

uo 

UO 

tin 

II 

1 u 

O 

0 

C 

0 

• 

M3 

r4 

• 

in 

CO 

•P 

IW 

CO  i 

• 

— 

• 

\ 

4-1 

r\ 

rjQ 

Mt 

O ! 

uo 

uo 

* 

P 

U 

CJ 

- 

- 

\ 

M3 

•r4 

1 

X 

0 

• 

CM 

O 

CO 

ON 

M3 

ON 

P'- 

M3 

co 

00 

t-4 

UO 

CO 

00 

x — 1 

Mt 

P^ 

P 

•4 

1 

O 

00 

M3 

CM 

ON 

Mt 

ON 

r-4 

O 

r-4 

GO 

00 

t— — 1 

CM 

p^ 

p"- 

M3 

CO 

4-1 

4-4 

r-4 

• 

PQ 

4-4 

I 

cd 

in 

M3 

00 

00 

p- 

in 

00 

00 

ON 

UO 

UO 

p-- 

00 

Mt 

ON 

r^. 

00 

CJ 

I ^ 

CJ 

in 

in 

uo 

UO 

uo 

uo 

uo 

UO 

uo 

UO 

UO 

in 

uo 

in 

UO 

uo 

in 

uo 

uo 

O 

i O 

p- 

O 

! 0 

• 

u 

• 

0 

00 

O 

U 

\ ■*-» 

CO 

• 

• 

CJ 

-a 

M3 

— 

-- 

-- 

-> 

: ON 

> ! 

II 

4-4 

O 

UO 

i in 

03 

j 

P 

I 

1 

1-4 

• 

CO 

00 

Mt 

uo 

ON 

M3 

Mt 

ON 

CM 

ON 

r-4 

r4 

ON 

r-4 

co 

Mt 

>N 

H 

| 

a 

CM 

CM 

I-'- 

Mt 

O 

uo 

ON 

M3 

i M3 

M> 

CM 

Mt 

CM 

Mt 

M3 

CO 

M3 

■Ul 

r-4 

•H 

4-1 

cd 

Mt 

UO 

UO 

M3 

UO 

Mi- 

M3 

uo 

M3 

UO 

M3 

r^. 

P" 

UO 

00 

p'- 

> 

cd 

! 

CJ 

in 

UO 

UO 

in 

uo 

UO 

UO 

uo 

in 

! uo 

in 

in 

uo 

UO 

uo 

uo 

UO 

in 

tM 

CJ 

T— 1 

4-1 

P3 

r— 1 

O 

cd 

P 

CJ 

IS 

• 

0 

O 

03 

a 

4-J 

CO 

• 



• 

i 

C 

cd 

r£3 

UO 

~ 

s 

ON 

O 

4-4 

i 0 

1 

uo 

UO 

CJ 

P 

P 

1 

r- 

0 

r^. 

0 

r-4 

CTJ 

CO 

uo 

O 

O 

in 

uo 

UO 

in 

uo 

CO 

i 10 

0 

0 

in 

UO 

uo 

UO 

uo 

co 

e 

M3 

cu 

fc-l 

u 

4-1 

0 <U 

i M3 

CM 

CM 

r-4 

0 

r-4 

CM 

CM 

O 

0 

CM 

CM 

r-4 

0 

r-4 

CM 

CM 

0 

0 

cj 

P 

OJ  0 

M3 

cd 

4-1  C 

H 

W 

U cd 

co 

co 

P 4J 

O 

O 

uo 

uo 

O 

O 

0 

Mf 

uo 

O 

O 

in 

UO 

0 

O 

O 

Mt 

UO 

CO  0 

Pd 

P 

M3 

r-4 

r-4 

r-4 

T— 1 

r-4 

r-4 

r-4 

r-4 

r-4 

r-4 

r-4 

r-4 

r-4 

r-4 

i-4 

T 1 

r-4 

r-4 

! 

P 03 

1 

cj  a 

p 0 

4 

1 — l 

r — 1 

1 

c 0 

T 

0 

0 

0 

O 

0 

in 

UO 

CM 

0 

0 

0 

0 

0 

0 

in 

UO 

CM 

0 

i 

M 

! £ 1 

| 

CM 

CM 

r- 4 

\ — 1 

r-4 

CM 

CM 

r-4 

r-4 

CM 

CM 

t-t 

r-4 

r-4 

CM 

CM 

r-4 

x — 1 

i 

II 

CJ 

o 

•H 


id 


0) 

pcj 


03 

03 

>0 

cd 

r-4 

cd 

03 

03 

Lewis 


NAPIER  SHELTER 


-35  - 


c\j 

od 

o u 
00  0 
C\J 

00 


C D 


3 

•H 


00 


CO 

o 


CD 

o 


St 

o 


00 

d 


•rH 


S3svo  jo  yjgiMnN 


INSIDE  SURFACE  HEAT  TRANSFER  COEFFICIENT  Btu/hrft 


- 3C  - 


Jo  ‘ JIAI31  JOVJdnSdJNNI  J3133HS  39VJ3AV 
- J1AI31  dIV  J3133HS  39Vd3AV  = IV 


SURFACE  HEAT  TRANSFER  COEFFICIENT,  Btu/hr  ft 


surface  heat  transfer  coefficients  upon  calculations  of  the  shelter 
thermal  environment  can  be  seen  in  table  2,  which  illustrates  a result 
of  sensitivity  analyses  of  inner  surface  heat  transfer  coefficients  using 
computer  program  M-(2)  on  the  NBS  shelter. 

In  table  2,  the  calculated  shelter  surface  temperatures,  shelter 
air  temperatures,  and  shelter  air  relative  humidities  are  compared  with 
the  observed  data  at  the  end  of  the  seventh  and  fourteenth  day  on  test 
5 of  the  NBS  shelter.  The  symbol  WG  in  table  2 is  zero  when  solar  heat 
input  to  the  ground  surface  above  the  shelter  is  ignored,  while  it  is 
one  when  the  entire  solar  heat  effect  is  assumed  absorbed  by  the  ground 
surface.  If  they  are  compared  at  an  identical  condition,  the  calculated 
shelter  air  temperature  at  the  fourteenth  day  became  approximately  3 F 
higher  when  solar  heat  input  was  considered  than  when  it  was  ignored. 

For  the  fourteenth  day  results,  the  interior  surface  temperature 
agreement  between  the  observed  and  calculated  is  better  with  rather  high 
values  of  surface  conductances,  while  better  air  temperature  agreement  is 
obtained  with  relatively  low  values  of  surface  conductance.  The  combina- 
tion of  h = 1.21,  h = 1.4^,  and  h = 0.57  is  obtained  from  reference 
w ’ R F 

L 10 j on  the  basis  that  the  simulated  occupants  obstruct  each  wall  from 

seeing  each  other.  As  far  as  the  shelter  air  temperature  is  concerned, 

a combination  of  h = 1.0,  ti  = 1.5,  and  h - 0.5  yielded  the  best 
w R ’ F 

agreement,  which  was  also  true  for  NBS  tests  3 and  4. 

5.6.  Outdoor  conditions. 

Outdoor  air  conditions  during  the  prototype  shelter  test  periods 

were  not  the  same  as  the  ventilation  air  conditions,  as  would  usually 

be  the  case  for  actual  non-air-conditioned  shelters.  The  ventilation  air 

conditions  were  selected  and  programmed  according  to  certain  climactic 

criteria  to  simulate  typical  operating  conditions  of  a shelter,  regardles 
of  the  actual  climatic  condition  during  the  test  period. 


- 38  - 


The  outdoor  air  temperatures  used  for  the  calculations  in  this 
report  were  recorded  separately  during  the  test,  averages  of  which  are 
shown  in  table  1. 

The  NBS  tests  included  observations  of  solar  energy  at  the  shelter 
roof  surface,  but  the  tests  conducted  by  the  University  of  Florida  did 
not  have  these  data.  Therefore,  the  inclusion  of  solar  energy  data  for 
the  shelters  tested  by  the  University  of  Florida  was  accomplished  by 
using  data  supplied  by  Flanigan  [ll]. 

5.7.  Initial  earth  temperature. 

As  described  in  section  3,  the  computer  model  M-(l)  incorporated 
an  earth  temperature  distribution  normal  to  the  surface  in  the  six 
surrounding  blocks.  The  NBS  test  shelter  included  these  observations  and 
they  were  used  for  the  M-(l)  calculation  of  this  shelter.  The  M-(2) 
calculation  provided  for  only  a vertical  or  depthwise  distribution  of 
initial  earth  temperature.  The  observed  data  on  roof  region  temperatures, 
average  of  the  wall  regions,  and  observed  distribution  of  floor  region 
were  used  to  arrive  at  a depthwise  distribution  of  the  earth  temperature. 
Average  earth  temperatures  of  the  wall,  roof,  and  floor  regions  were  used 
for  the  calculations  employing  M-(3)  and  M-(4)  models.  Table  1 shows  the 
overall  average  values  of  initial  earth  temperature  which  were  used  for 
the  calculations  in  these  latter  two  models*. 

5.8..  Comparison  between  the  computed  and  observed  shelter 
thermal  environments 

Figures  5 through  25  represent  some  of  the  results  obtained  by 
computer  analysis  of  thermal  environment,  together  with  the  observed 
data.  Figures  5,  6,  7,  and  8 compare  the  calculated  shelter  thermal 
environments  with  the  observed  data  for  tfcsts  1,  2,  and  3 on  the  NBS 
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program  M-(l)) . 
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Fig.  8 Comparison  of  the  calculated  and  observed  shelter  air  temperatures 
and  relative  humidity  for  test  3 of  NBS  family  shelter  (computer 
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family  shelter.  The  M-(l)  computer  program  was  used  for  these  compari- 
sons. Computed  shelter  air  temperatures  agreed  with  the  observed  within 
approximately  2 °F  for  the  entire  period  of  shelter  occupancy.  Agreement 
between  computed  and  observed  relative  humidities  was  good  for  tests  1 
and  3,  and  somewhat  poorer  for  test  2 using  the  M-(l)  program.  Figure  7 
represents  a sealed-up  condition  with  6 occupants,  as  indicated  in  table  1. 
A very  good  agreement  in  calculated  and  observed  air  temperature  for  the 
first  4-day  period  is  shown  in  figure  7.  However,  the  calculated  relative 
humidity  was  as  much  as  6 percent  lower  than  the  observed  during  the 
second  day  of  test,  and  did  not  quite  attain  the  saturated  condition  that 
was  observed  during  the  test,  even  after  48  hours  of  sealed-up  condition. 

Figures  9,  10,  and  11  compare  temperatures  and  humidities  obtained 
by  computer  program  M-(2)  with  the  observed  conditions  for  NBS  shelter 
tests  3,  4,  and  5.  Figures  8 and  9,  both  for  test  3 condition,  show  that 
computer  program  M-(l),  the  temperature-block  model,  results  in  better 
agreement  with  the  observed  shelter  temperatures  for  the  first  two  days 
than  that  by  program  M-(2),  the  homogeneous -earth  model.  This  was 
expected,  because  program  M-(l),  as  explained  before,  is  considerably 
more  elaborate  in  accounting  for  a complex  nature  of  initial  earth  tem- 
perature distribution  around  the  shelter  than  program  M-(2),  and  also 
takes  into  account  the  differences  in  thermal  properties  of  the  concrete 
and  the  earth. 

A.  However,  the  agreement  between  the  calculated  and  observed  sir 
temperature  at  the  5-foot  level  in  the  shelter  was  somewhat  better  using 
the  M-(2)  program  than  the  M-(i)  program  during  the  last  4 or  5 days  of 


test  3. 


- 44  - 


The  agreement  of  the  computed  shelter  air  temperature  with  the 
observed  condition  for  test  4 of  the  NBS  shelter  became  poorer  toward  the 
end  of  the  test.  The  lowering  of  the  computed  shelter  temperature  during 
the  second  week  of  the  test  was  caused  primarily  by  a decrease  in  outdoor 

temperature.  This  temperature  decrease  averaged  approximately  15  degrees, 
beginning  on  October  12  and  continuing  to  the  end  of  the  test  period.  The 
detailed  report  [3]  on  these  tests  indicates  a steady  and  significant 
decline  of  earth  temperature  during  this  period.  The  curves  on  figure  10 
shows  that  the  observed  shelter  air  temperature  was  not  affected  as  sig- 
nificantly by  this  cool  spell  as  the  computer  model  indicated  that  it 

wou Id  be  . 

The  sudden  and  irregular  drop  of  the  calculated  relative  humidity 
for  the  NBS  test  5 condition  shown  in  figure  11  is  not  reflected  by  the 
observed  data.  It  is  probable  that  the  relative  humidity  in  the  shelter 
was  sustained  at  a high  level  by  drying  of  the  shelter  walls.  The  mois- 
ture balance  between  the  supply  and  exhaust  air  indicated  that  this 
evaporation  amounted  to  approximately  2 to  2 .5  lb/day  during  this  winter 
test  condition  where  extremely  dry  ventilation  air  was  employed.  As 
indicated  earlier,  the  computer  program  developed  for  this  analysis  had 
no  provision  for  taking  this  drying  process  into  consideration.  Table  2 
indicates  that  the  calculated  air  temperature  at  the  end  of  14  days  in 
the  NBS  test  shelter  test  5 condition  is  in  better  agreement  with  the 
observed  value  when  solar  heat  effect  was  not  included  for  the  ground 
surface  heat  exchange  than  when  it  was  included.  The  ground  surface  was 
quite  wet  during  this  test  period  because  of  a snow  prior  to  the  test. 

This  implies  that  during  this  particular  test  period  the  solar  heat  was 
mostly  absorbed  by  evaporation  of  water  on  or  near  the  earth* s surface 
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Fig.  9 Same  as  Fig.  8 but  using  computer  program  M- (2) . 
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relative  humidities  for  test  4 of  NBS  family  shelter  (computer 


47 


05 

u 

<0 

<o 

U 

•u 

3 

u 

a 

CO 

£ 

Pi 

o 

CO 

o 

& 

'w' 

0) 

Pi 

OJ 

co 

•u 

Pi 

i — i 

<0 

CO 

05 

Pi 

CO 

>> 

4J 

i — 1 

•rl 

0) 

& 

cO 

05 

4-1 

to 

C/0 

co 

PQ 

> 

S3 

u 

(1) 

4-1 

05 

O 

O 

in 

•tJ 

u 

£ 

CO 

CO 

CO 

■U 

t) 

0) 

u 

4J 

o 

CO 

4-! 

H 

2 

05 

O 

<0 

i — 1 

-H 

cO 

4J 

o 

TO 

<0 

°r^ 

& 

g 

U 

• 

<4-1 

o 

CO 

<S 

> 

v_^ 

a 

t 

o 

4J 

S 

05 

cO 

•H 

1 — I 

g 

co 

cO 

CO 

p< 

Pi 

a. 

50 

& 

rO 

O 

o 

pi 

Pi 

o 

cO 

a. 

r-t 

1 t 

» 

50 

•H 

P4 

- 48  - 

and  did  not  cause  appreciable  temperature  rise  in  the  ground.  Similar 
moisture  evaporation  near  the  earth's  surface  also  explains  the  small 
earth  temperature  rise  near  the  surface  for  other  test  conditions,  such 
as  the  Summerlin,  Broyles,  and  Napier  shelters,  some  of  which  will  be 
illustrated  later. 

Figure  12  compares  calculated  (M-(2)  program)  and  observed  air 
temperatures  in  the  Summerlin  shelter  for  its  summer  condition  test. 
Although  the  overall  trend  of  the  actual  temperature  is  closely  followed 
by  the  calculation,  the  computed  amplitude  of  the  diurnal  temperature 
variation  within  the  shelter  is  considerably  smaller  than  the  observed. 
The  relative  humidity  calculation  for  the  summer  test  condition  of  the 
Summerlin  shelter  is  shown  in  figure  13.  For  the  first  five  days,  the 
calculated  relative  humidity  was  higher  than  the  observed.  This  dis- 
crepancy may  be  due  to  the  fact  that  there  were  cool  regions  in  the 
shelter  interior  surfaces  where  more  condensation  of  water  vapor  was 
taking  place  than  was  determined  by  the  calculation  which  was  based  upon 
average  surface  temperatures  of  each  exposure . Some  of  the  high  peaks  in 
the  computed  relative  humidity  during  the  last  three  days  were  not  regis- 
tered in  the  observed  record.  Figure  14  compares  the  calculated  earth 
temperature  surrounding  the  Summerlin  shelter  during  the  summer  test 
period  with  the  observed.  The  agreement  between  calculated  and  observed 
values  was  quite  good  for  the  south  and  east  walls.  The  earth  temperature 
under  the  floor  was  not  observed  during  the  test  because  this  was  a pri- 
vately owned  shelter  and  it  was  decided  that  the  water-tightness  of  the 
floor  should  not  be  jeopardized  by  making  a hole  through  it.  The  cal- 
culated roof  region  temperature  was  much  higher  than  the  observed  values. 
As  mentioned  before,  this  discrepancy  is  assumed  to  be  caused  by  partial 
utilization  of  solar  energy  by  surface  evaporation  of  ground  moisture, 
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Fig.  13.  Comparison  of  the  calculated  and  observed  shelter  air  relative  humidities 
for  summer  condition  test  for  Summerlin  Shelter  (computer  program  M- (2)  ). 
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of  the  calculated  and  observed  earth  temperatures  surrounding 
helter  during  summer  test  conditions  (computer  program  M-(2) . 
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resulting  in  less  temperature  rise  in  the  interior  of  the  earth  than  cal- 
culated. It  was  found  for  all  other  shelters  for  which  the  ground  surfac 
was  not  dry  and  bare,  the  inclusion  of  solar  heat  in  the  calculation 
usually  gave  higher  shelter  temperatures,  as  well  as  higher  roof  region 

temperatures,  than  observed  values. 

Figure  15  compares  the  computed  Summerlin  shelter  thermal  environ- 
ment with  the  observed  condition  for  a moderate  climatic  condition. 
Although  the  computed  temperature  amplitude  was  again  lower  than  the 
observed,  the  agreement  between  the  computed  and  the  observed  data  is 
very  good  and  was  better  than  for  the  hot  summer  test  condition.  Three 
high  peaks  of  the  observed  shelter  relative  humidity  during  the  last  two 
days  of  the  test  were  not  reproduced  by  the  calculation  and  they  may 
represent  instrumentation  errors . Essentially  identical  results  were 
obtained  when  the  calculations  were  repeated  by  M-(3)  or  one -dimensional 
compound  model  on  Summerlin  shelter. 

Figures  16  and  17  compare  calculated  results  with  observed  tem- 
perature and  relative  humidities  of  the  Broyles  shelter.  As  described  in 
section  4-C,  this  shelter  was  conditioned  by  a cooling  coil  using  well 
water  during  the  first  two  days.  In  figure  16,  a dashed  curve  shows  the 
computed  result  without  consideration  of  the  cooling  coil,  thus  yielding 
Tauch  higher  shelter  temperatures  during  the  first  two  days  of  the  test. 
The  last  2-day  portion  of  figure  16  was  not  simulated  by  the  computer, 
so  the  comparison  between  the  calculated  and  the  observed  shelter  tem- 
perature without  air  conditioning  should  be  made  between  August  1 and  16. 
However,  an  adjustment  was  made  later  to  account  for  the  observed  cooling 
capacity  of  the  coil  by  making  Q^S  an<^  Qml  ne8ative  the  equation  for 

miscellaneous  heat  load  in  section  2.  The  first  2-day  portion  of  the 
Broyles  shelter  calculation  was  repeated  with  this  adjustment  and  the 
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results  are  shown  by  solid  dots  on  figure  16.  The  adjusted  and  cal- 
culated shelter  air  temperatures  for  this  initial  2-day  period  agree 
very  well  with  observed  temperatures  obtained  at  the  geometric  center  of 
the  shelter 

Figure  17  compares  the  calculated  and  observed  relative  humidi- 
ties in  the  Broyles  shelter.  The  large  disagreement  of  the  computed 
relative  humidities  during  August  2 to  7 cannot  be  explained  adequately 
by  the  available  data.  The  agreement  of  the  computed  relative  humidities 
with  the  observed  values  for  this  test  were  generally  very  poor,  and  the 
same  trend  was  also  observed  in  the  Napier  shelter  comparison,  as  shown 
in  figure  19. 

The  Napier  shelter  temperature  comparison  was  satisfactory  if  the 
solar  heat  input  to  the  ground  surface  is  ignored,  as  seen  from  figure  18. 
The  inclusion  of  the  solar  heat  in  the  calculation  caused  the  computed 
shelter  air  temperature  to  be  approximately  3 to  4 degrees  higher  than  the 
observed  temperature . 

By  contrast,  figure  20  shows  excellent  agreement  between  the 
calculated  and  observed  temperature  and  relative  humidity  during  the 
winter  test  condition  for  the  Reading  shelter.  For  this  shelter,  nearly 
80  percent  of  the  total  heat  generated  was  conducted  into  the  surrounding 
earth.  The  moisture  balance  of  the  test  results  is  shown  in  figure  21, 
indicating  the  continuous  condensation  of  water  vapor  on  the  inner  sur- 
faces of  the  shelter  during  the  first  ten  days  of  the  test.  The  venti- 
lation rate  was  either  3 cfm  or  \\  cfm  per  person  during  this  period. 

The  excellent  agreement  between  the  computed  thermal  environment  with 
the  observed  condition  during  the  first  10-day  period  is  a good  demon- 
stration that  the  computer  simulation  technique  was  valid.  However,  in 
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figure  20,  the  computed  relative  humidity  is  approximately  10  percent 
lower  than  the  observed  during  the  last  three  days.  During  this  period, 
the  simulated  occupants  were  increased  from  the  equivalent  of  50  to  100. 

In  figure  21,  the  indicated  moisture  balance  within  the  Reading  shelter 
during  this  period  shows  that  the  ventilation  air  actually  carried  out 
more  water  vapor  than  was  produced  by  all  of  the  simulated  occupants. 

This  proves  that  the  shelter  wall  was  drying  out  during  the  period.  As 
mentioned  before,  the  computer  program  does  not  have  provision  for  simu- 
lating the  drying  process  of  concrete  walls;  therefore,  the  computed 
relative  humidity  was  lower  than  the  observed. 

Figures  22  and  23  compare  the  computed  and  observed  earth  tempera- 
tures surrounding  the  Reading  shelter.  Being  in  the  midst  of  winter, 
the  undisturbed  earth  temperature  was  relatively  low  and  the  ground  sur- 
face was  snow  covered  during  this  period.  Two  computer  models  M-(2) 
(symmetrical  and  three-dimensional)  and  M-(3)  (one -dimensional  and  compound), 
were  applied  to  the  Reading  shelter  calculation,  resulting  in  practically 
identical  thermal  environments.  The  ground  temperature  profiles  were 
computed  and  compared  with  the  observed  profiles  for  other  prototype 
shelters  and  the  agreement  between  the  calculated  and  observed  results 
were  generally  similar  to  those  shown  in  figures  22  and  23  for  the 
Reading  shelter. 

The  computer  Model  M-(4),  basically  a one -dimensional  and  compound 
system  with  the  top  surface  of  the  roof  being  adiabatic,  was  employed  to 
compute  the  thermal  environments  of  the  Ft.  Belvior  200-man  and  1000-man 
basement-type  shelters.  These  results  are  shown  graphically  in  Figures  24 
and  25.  Instead  of  comparing  the  relative  humdity,  as  in  cases  of  other 
shelters,  effective  temperatures  were  compared,  in  addition  to  the  dry-bulb 
temperature  comparison.  The  computed  and  observed  results  agreed  almost 
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Fig.  21  Observed  moisture  balance  of  Reading  Shelter. 
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Fig.  22  Comparison  of  the  calculated  and  observed  earth  temperatures  outside 
Reading  Shelter  walls  (computer  program  M-(2) ) . 
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HEADING  SHELTER  EARTH  TEMPERATURE  ON 
Nth  DAY  OF  THE  TEST 
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Comparison  of  the  calculated  and  observed  earth  temperature  outside 
Reading  Shelter  ceiling  and  floor  (computer  program  M- (2) ) . 
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CALENDER  DATE,  SEPTEMBER  AND  OCTOBER,  1963 

Comparison  of  the  calculated  and  observed  shelter  air  temperatures  and 
effective  temperatures  of  Ft.  Belvoir  200-man  shelter . (computer  program 
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perfectly  for  the  1000-man  shelter.  Higher  and  less  fluctuating  dry- 
bulb  and  effective  temperatures  were  calculated  for  the  200-man  shelter 
than  were  observed,  particularly  during  the  first  seven  days  of  the  test 
period.  By  comparison  it  appeared  that  the  observed  temperatures 
reflected  shelter  diurnal  cycles  corresponding  to  a higher  ventila- 
tion air  rate  than  the  indicated  3 cfm  Per  person.  It  also  should  be 
realized  that  the  geometrical  relationship  between  the  inlet  for  the 
ventilation  air  and  the  stations  for  room  temperature  measurement  and 
the  mixing  effects  in  the  intervening  space  would  have  an  effect  in  the 
amplitude  of  the  temperature  variations  observed  in  any  of  the  shelters. 

6.  CONCLUSIONS 

Except  for  a very  few  cases,  such  as  for  the  Napier  shelter  and 
the  NBS  shelter  test  4 conditions,  the  computer  simulation  based  upon  a 
finite  difference  solution  of  the  heat  conduction  equation  for  the  cal- 
culation of  the  thermal  environment  of  underground  protective  shelters 
has  been  found  generally  satisfactory.  Good  agreement  between  calculat- 
ed and  observed  thermal  conditions  was  obtained  for  tests  2 and  3 i-n  the 
NBS  shelter,  the  Summerlin  shelter,  the  Reading  shelter,  and  the  Ft. 
Belvoir  1000-man  shelter.  This  agreement  coincided  with  well-controlled 
test  conditions,  and  where  input  parameters  were  known  more  accurately 
than  in  other  cases.  However,  for  this  type  of  analysis  there  are  two 
inherent  uncertainties  which  influence  the  accuracy  of  the  calculations. 

The  first  uncertainty  stems  from  simplification  of  the  actual 
environmental  system  when  constructing  the  mathematical  model.  The 
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complexity  of  the  system  including  construction  details  of  interior 
partitions,  foundations,  wall  supports  and  reinforcing;  the  heter- 
ogeneous temperatures,  thermal  properties,  and  types  of  the  surround- 
ing earth;  and  space  variation  of  conditioning  air  and  its  convective 
pattern  in  the  shelter,  makes  an  accurate  mathematical  description 
of  the  system  very  difficult.  Therefore,  the  simulated  computer 
model  is  a simplified  version  of  the  actual  system  that  is  suited  to  the 
mathematical  handling  of  the  problem.  The  magnitude  of  the  error  due  to 
this  uncertainty  has  not  been  evaluated,  but  it  is  a part  of  the  differ- 
ences between  calculated  and  observed  temperature  and  humidity  conditions 
shown  in  this  report  for  seven  different  shelters. 

The  second  uncertainty  is  closely  related  to  the  first,  and  con- 
cerns the  reliability  of  input  data.  It  is  almost  impossible  to  acquire 
complete  three-dimensional  information  on  earth  thermal  properties  and 
temperature,  and  if  such  information  were  available,  it  would  be  difficult 
to  use  as  computer  imput  data.  Therefore,  when  a homogeneous  heat  con- 
duction model  is  employed  for  a system  surrounded  by  earth  of  hetero- 
geneous characteristics  and  temperature,  the  choice  of  proper  input  values 
becomes  very  difficult.  The  uncertainty  involved  in  evaluating  the  interior 
surface  heat  transfer  coefficients  is  equally  as  difficult,  as  discussed 
in  section  5. 

Thus  the  simplification  of  the  physical  characteristics  of  the  real 
shelter-earth  system  employed  in  creating  an  analytical  model  and  the 
uncertainties  in  the  thermal  properties  of  the  materials  involved, 
became  the  principal  approximations  incorporated  into  the  computer 
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analysis  of  the  thermal  environment  in  shelters.  In  spite  of  these 
uncertainties , good  agreement  between  computed  and  observed  results  was 
achieved  for  many  of  the  shelters  descirbed  in  this  report. 

In  addition  to  these  general  comments,  these  specific  conclusions 
are  emphasized,  as  follows: 

(a)  Soil  analysis  of  most  of  the  prototype  shelters  indicat- 
ed that  thermal  diffusivity  and  thermal  conductivity  are  in  the  neighbor- 
hood of  0.02  ft8 /hr  and  0.75  Btu/hr  ft2 (deg  F/ft) , respectively.  These  values 
in  turn,  seem  to  result  in  fairly  good  agreement  of  earth  temperature 

change  during  the  test  period  for  most  of  the  shelters. 

(b)  Heat  conductance  at  the  shelter  inner  surfaces  for  most 

f 

of  the  summer  shelter  conditions  may  be  approximated  by  the  following: 

1.0  Btu/hr  ft2  deg  F for  vertical  walls. 

1.5  Btu/hr  ft2  deg  F for  the  ceiling. 

0.5  Btu/hr  ft2  deg  F for  the  floor. 

Although  there  may  be  some  other  values  and  other  combinations  of  these 
values  that  may  have  resulted  in  a slightly  better  simulation  than  those 
used  in  the  analysis,  these  three  values  can  be  considered  representative 
design  heat  transfer  coefficients  in  the  underground  cavities. 

(c)  For  large  shelters,  the  one-dimensional,  compound  6- 
directional  model  M-(3)  will  probably  be  adequate  for  calculating  the 
shelter  heat  transfer.  Thus,  the  complicated  3 -dimensional  model  may  not  be 
required,  at  least  for  the  heat  transfer  calculation  of  less  than  14-day 

\ J 

occupancy.  For  small  shelters  (such  as  a family  shelter  similar  to  the  NBS 

shelter) , however,  it  is  recommended  that  the  three-dimensional  model  be 

used  for  the  accurate  calculations. 
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7 . APPENDIX 

7 . 1 Air-conditioning  subroutine. 

Program  M-(4)  employs  a subroutine  for  air  conditioning  the 
shelter.  This  subroutine  enables  the  computer  program  to  include  the 
heat  and  vapor  absorbing  capacity  of  a given  cooling  system  (a  fan-and- 
coil  unit)  in  the  shelter  heat  balance  equation  of  2.1.  The  row  by  row 
account  of  the  cooling  coil  performance  is  considered  for  the  cross- 
counter flow  circuited  coil  circulating  the  well  water.  The  input  data 
required  for  this  subroutine  are: 

Shelter  air  dry-bulb  temperature,  t °F 

Shelter  air  wet-bulb  temperature,  t^,  °F 

Ventilation  air  dry-bulb  temperature,  t 


v 

Ventilation  air  wet-bulb  temperature,  T , 


Inlet  coolant  temperature,  t 


w 


Recirculation  air  rate,  CFM^ 

Contact  factor  of  air  conditioning  coil  per  row,  C 
Thermal  resistance  between  the  solid-air  inter- 


face and  the  coolant,  R 


w 


Estimated  temperature  rise  due  to  fan  heat,  At 


Coolant  heat  content,  G 


w 


F,  ft2  hr/Btu 


Btu/hr,  °F 


The  contact  factor  per  row  C^,  mentioned  above,  is  a function  of  air 


face  velocity  across  the  given  air  cooling  coil,  air-side  heat  transfer 
coefficient  of  the  coil  surface,  and  the  amount  of  total  heat  transfer 
surface.  The  factor  may  be  estimated  by  the  following  expressions: 
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where f = air  side  heat  transfer  coefficient,  Btu/hr  ft~  °F. 
a 

S = air  side  heat  transfer  surface  of  the  coil  per  row,  ft 
P 

CFMt  = CFM^  + CFM  = total  air  flow  rate. 

The  thermal  resistance  value  R may  be  estimated  by  the  following  expression 


JL 

R 

w 


S 

_£ 


+ Hr 


7-2 


where  St  = coolant  side  heat  transfer  surface  of  the  coil  per  row,  ft2 
f = heat  transfer  coefficient  of  coolant,  Btu/hr,  ft2  °F. 

Hr  = sum  of  all  other  heat  resistance  between  the  coolant  side 
surface  and  air  side  surface,  such  as  thermal  resistance 
due  to  tube  wall,  finish  bond,  finish  metal  (depending  upon 
the  effectiveness  of  finish),  and  condensate  film  thickness. 

The  temperature  rise  due  to  fan  heat  can  be  estimated  by  total 
power  input  to  the  fan  and  air  flow  rate,  assuming  that  all  the  fan  heat 
will  be  expended  to  raise  the  air  stream  temperature.  The  value  of  Atf 
usually  does  not  exceed  3 °F. 

The  coolant  heat  Gw  is  the  coolant  mass  flow  rate  multiplied  by 
the  coolant  specific  heat.  In  the  case  of  a direct  expansion  refrigerant 
coil,  where  the  coolant  temperature  change  in  the  coil  is  very  small,  Gw 
is  considered  infinity,  or  a very  large  number.  The  details  of  the  cal- 


71 


culative  procedure  on  Cf  and  Rw,  and  their  design  value  can  be  found  in 
references  5 and  6. 

First,  the  air  inlet  condition  to  the  air  cooling  coil  is  calculated 
by  the  following  equation: 

(CFM)R(ta)  + (CFM)(tv) 

t±  ~ CFMr  + CFM  7-3 

_ (CFM)R(Wa)  + (CFM)(WV) 

1 “ (CFM)r  + CFM  7-4 

where  t-j_,  Wi?  Wa,  and  Wy  represent  inlet  air  dry-bulb  temperatures,  inlet 
air  humidity  ratio,  shelter  air  humidity  ratio,  and  ventilation  air 
humidity  ratio,  respectively.  These  humidity  ratio  values  are  calculated 

by  a separate  subroutine  from  dry-  and  wet-bulb  temperatures,  as 
explained  in  section  2.7. 

The  air  conditioning  coil  is  assumed  to  be  of  multi-row  structure, 
and  the  overall  direction  of  coolant  is  counter  to  that  of  air  flow, 
although  for  individual  rows  the  coolant  is  flowing  perpendicular  to  the 
air  stream.  For  simplicity,  it  is  assumed  that  the  coolant  temperature, 
as  well  as  air-side  surface  temperature  of  a row,  changes  by  a step  func- 
tion. The  heat  and  vapor  transfer  calculation  of  the  counter-flow  dehumidi- 
fying  coil  requires  an  iterative  procedure,  because,  except  for  the  direct 
expansion  coil,  the  outlet  condition  of  air  or  coolant  is  not  previously 
known.  In  this  report,  the  outlet  coolant  temperature,  which  is  in  the 
same  side  as  the  inlet  air  condition  of  the  coil,  is  first  approximated. 

The  heat  and  vapor  transfer,  and  subsequent  reduction  of  air  temperature, 
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air  humidity,  and  change  of  the  coolant  fluid,  are  then  calculated 
row  by  row.  After  the  calculation  is  completed,  the  total  net  change 
in  coolant  temperature  is  subtracted  from  the  outlet  temperature  of 
the  coolant,  yielding  the  calculated  coolant  temperature  at  the  inlet 
condition.  If  the  calculated  coolant  temperature  at  the  coil  coolant 
inlet  is  different  from  the  actual,  the  calculation  is  repeated  with  a 
modified  outlet  coolant  temperature  until  the  calculated  agrees  with  the 
given  inlet  temperature  of  the  coolant  flow.  The  general  heat  and  vapor 
transfer  relation  used  for  this  calculation  is  described  for  i row  as 
follows : 


WS£  in  the  above  equation  is  the  humidity  ratio  of  the  air  saturated  at 
the  surface  temperature,  ts-j_;  if  ^ Ws^,  t^ie  second  term  in  the  left 
hand  side  of  the  equation  is  set  equal  to  zero.  Since  Ws^  is  a compli- 
cated function  of  tg  ^ , an  iterative  technique  is  required  to  solve  ts^ 
from  the  above  expression.  After  ts^  is  obtained,  the  leaving  air  condi- 
tion, tj[+q,  W^+]_,  from  the  i row,  and  entering  coolant  condition  to  the 
i + 1,  are  calculated  by  the  following  relations: 


7-5 


ti  “ ti+l  (fci  “ tsHCf) 


7-6 


Wi  - Wi+l  — (Wj_  - Wgi)(Cf)  if  W£  > Wg£ 


7-7 


0 


— ^-(tsi  - twj_)  - (twj_  - tWj_-(-i)(Gw) 


7-8 


Since  all  of  the  relations  are  linear,  calculations  of  t^+^,  and 

tWji_i  are  straightforward. 

When  the  total  number  of  rows  is  N and  the  properties  in  the  air 
inlet  point  are  specified  by  subscript  1,  the  properties  of  the  air  out- 
let condition  should  be  specified  by  N + 1.  Thus,  the  calculations  will 
be  iterated,  as  mentioned  before,  until  tW}^+l  becomes  equal  to  tw . The 
final  results  of  the  air  conditioning  capacity  will  be  expressed  as 
follows : 

(1)  Sensible  cooling  capacity,  Q^g  = (1 .08)  (CFM)  (t^-t^^+Atf ) 

(2)  Latent  cooling  capacity,  Qal  = (4  .5  ) (CFM)  (W]_-Wjj+]_)\  . 

With  the  value  of  Q^g  and  Q^l  known,  the  shelter  air  condition 

for  the  next  time  period  now  can  be  calculated  by  adding  Q^s  and  Q^l  to 
the  overall  heat  balance  equation  in  (1). 
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Computer  Programs 


Fortran  listings  of  the  computer  programs  developed  for  the 
heat  transfer  studies  in  underground  protective  shelters  during  the 
contract  research  are  attached  in  the  following  pages.  Except  for 
M-3  program,  the  computer  symbols  are  explained  according  to  input 
and  output  sequences  at  the  beginning  of  each  program. 

A complete  set  of  input  data  and  illustrative  examples  of 
output  format  for  program  M-4  are  also  shown  . 
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Program  M-l  Input  and  Output  Symbols 


KPl: 

A matrix 

index  marking 

the 

interface  of 

earth 

KP2: 

A matrix 

index  marking 

the 

interface  of 

earth 

KP3: 

A matrix 

index  marking 

the 

end  of  earth 

block 

wall . 

KQl: 

A matrix 

index  marking 

the 

interface  of 

ear  th 

KQ2: 

A matrix 

index  marking 

the 

interface  of 

earth 

KQ3: 

A matrix 

index  marking 

the 

end  of  earth 

block 

wall . 

KRl: 

A matrix 

index  marking 

the 

interface  of 

ear  th 

KR2 : 

A matrix 

index  marking 

the 

interface  of 

earth 

KR3: 

A matrix 

index  marking 

the 

end  of  earth 

block 

floor . 

KP 

KQ 

KR 

NN 

LC 

WC: 


TGI: 


Nearest  integer  to  (KPl  + KP2)/2. 

Nearest  integer  to  (KQl  + KQ2)/2. 

Nearest  integer  to  (KRl  + KR2)/2. 

Number  of  time  increments. 

Number  of  matrix  points  for  shelter  inner-wall. 

Work  cell  for  determining  type  of  initial  earth  temperature  data, 

WC  = -1.  : Depthwise  variation 

= 0.  : Homogeneous  earth  temperature 

=1.  : Six  directional  variations. 

A constant  initial  earth  temperature  when  WC  = 0. 


C-2 


TG1(I),  TG2(I),  TG3(I) , TG4(I),  TG5(I),  TG6(I):  Initial  earth  tempera- 

ture profiles  in  earth  blocks  outside  of  north,  south,  east, 


west,  roof  and  floor,  respectively.  °F 

TC(L,M) : 

Initial  temperature  distribution  of  the  shelter  inner  wall, 

S(M) : 

°F  o 

2 

Shelter  inner  wall  surface  area,  ft  . 

D(M): 

Thickness  of  shelter  inner  wall,  ft. 

ZL(M): 

Thickness  of  the  earth  block,  ft. 

CK(M): 

Thermal  conductivity  of  shelter  wall,  Btu/hr,ft,°F 

CG(M): 

Thermal  conductivity  of  earth  block,  Btu/hr , ft , °F 

AC  (M) : 

2 

Thermal  diffusivity  of  shelter  wall,  ft  /hr. 

AG(M): 

2 

Thermal  diffusivity  of  earth  block,  ft  /hr. 

H(M) : 

Heat  transfer  coefficient  at  the  shelter  inner  surface, 
Btu/hr , ft2 , °F. 

HD(M): 

Mass  transfer  coefficient  at  the  shelter  inner  surface, 
lb/hr,  ft2(lb/lb). 

HG(H) : 

Heat  transfer  coefficient  at  the  external  side  of  the  earth 

2 

block,  Btu/hr,  ft  , °F. 

DB(N) : 

Outdoor  air  dry-bulb  temperature,  °F. 

TV(CN) : 

Ventilation  air  dry-bulb  temperature,  °F. 

DPV(N) : 

Ventilation  air  dew-point  temperature,  °F. 

Q(N): 

2 

Solar  radiation*-  Btu/hr, ft  . 

A; 

Shelter  length  (internal  dimension),  ft. 

B: 

Shelter  width  (internal  dimension),  ft. 

C: 

Shelter  height  (internal  dimension),  ft. 

C-3 


ZN1: 

Number  of  600  Btu/hr  occupants. 

ZN2: 

Number  of  400  Btu/hr  occupants 

ZN3: 

Number  of  200  Btu/hr  occupants 

TA: 

Shelter  air  temperature,  °F. 

DPA: 

Shelter  air  dew-point  temperature,  °F. 

PB: 

Barometric  pressure,  in.  Hg . 

BS: 

Total  sensible  heat  emitted  by  non-human  heat  source,  Btu/hr 

BL: 

Total  latent  heat  generated  by  non- human  heat  source,  Btu/hr. 

DGX: 

Finite  Difference,  earth  length  for  N-S.  Coordinate  direction, 
ft. 

DGY: 

Finite  Difference,  earth  length  for  E-W.  Coordinate  direction, 
ft. 

DGZ: 

Finite  Difference,  earth  length  for  roof-floor  coordinate 
direction,  ft. 

DGX1: 

Modification  of  DGX  for  KPl  <1  <KP2,  ft. 

DGYl : 

Modification  of  DGY  for  KRl  <J  <KR2,  ft. 

DT: 

Finite  difference  inner  wall  time,  hr. 

V: 

Total  interior  volume  of  the  shelter,  cu  ft. 

G: 

Ventilation  air  rate,  6fm. 

ZN: 

Run  identification  number. 

DTT: 

Finite  difference  time  for  earth  block,  hr. 

TG(I,J,K) 

: Earth  temperature  °F. 

N: 

Number  of  time  iterations  (subscript). 

TM: 

Elapsed  time,  hr. 

QVS: 

Sensible  Heat  Exchanged  with  Ventilation  Air,  Btu/hr. 
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QVL: 


Latent  Heat  Exchanged  with  Ventilation  Air,  Btu/hr. 

QGS:  Sensible  Heat  Generated  within  the  shelter,  Btu/hr. 

QGL:  Latent  Heat  Generated  within  the  shelter,  Btu/hr 

QWST:  Heat  Transferred  to  the  shelter  inner  surfaces,  Btu/hr. 

QWLT:  Latent  Heat  Transferred  to  the  shelter  inner  surfaces,  Btu/hr. 

TWCT:  Water  Vapor  Collected  on  the  shelter,  inner  surfaces,  lb/hr. 

TQVS , TQVL,  TQGS,  TQGL , TQWST,  TQWCT , TTWCT:  Cumulative  values  for 

QVS,  QVL,  QGS,  QGL,  QWST,  QWLT  and  TWCT,  respectively. 

AT(N)  = TA  Calculated  shelter  air  dry -bulb  temperature,  °F. 

RH(N)  = RHA:  Calculated  shelter  air  relative  humidity,  %. 

DP(N)  = DPA:  Calculated  shelter  air  dew-point  temperature,  °F. 

Subscript  M refers  to  exposures  such  that 
M = 1:  North  wall  and  its  region 

2:  South  wall  and  its  region 

3:  East  wall  and  its  region 

4:  West  wall  and  its  region 

5:  Roof  and  its  region 

6:  Floor  and  its  region. 

Subroutine : 

GN : Computes  the  moisture  saturated  air  vapor  pressure  for  a given 

temperature  using  Goff  and  Gratch  formula. 

FN : Computes  the  boundary  temperature  by  a finite  difference  formula. 

QG : Computes  the  human  metabolic  heat  by  type  of  occupants  as  function 
of  temperature. 


C-5 


o o 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


UNDERGROUND  FALLOUT  SHELTER  THERMAL  ENVIRONMENT  HEAT  TRANSFER  1185000C 

NATIONAL  BUREAU  OF  STANDARDS  PRO JECT- 1 0436  T.KUSUDA,  10.3  I185001C 

DIMENSION  TG(23,22,29),TC( 10,6 ) ,QSUN( 1 70 , 6 ) , DB ( 170 ) , TV ( 170 ) 

DP V { 170  ) , TGH  10)  , TG2U0)  , TG31 10)  »TG4(10)  ,TG5  (10)  1185003C 

TG6(  10)  ,TG0(30)  , S ( 6 ) ,D  { 6')  , ZL  ( 6)  , AG  16  ) , AC  (6  ) , CK  (6  ) , CG ( 6 ) 1 1 85004C 
H(6) fHD(6) ,HG(6) ,QWS(6) ,QWL(6) , TQWS(6) , TQWL(6) ,TWC(6)  1185005C 

TTWC(6),QSUNT(6)  , AT( 170) » RH  ( 170) pDP( 170)  ,QSUN0(6) 

DIMENSION  CD(6) , E ( 6 ) , Q ( 170) 

3 FORMAT { 8F9. 4 ) 

1 FORMAT ( 14I4t2F3.0) 

99  READ  1,KP1,KP2»KP3,KQ1,K02,KQ3,KR1,KR2,KR3,KP,KQ,KR,NN,LC,HC 
IF(WC)  2,  12,11 
12  READ  3, TGI 
GO  TO  10 

2 GO  TO  13 

11  READ  3, (TGH I ) , I=1,KP1)» (TG2( I ) »I=KP2»KP3) , ( TG3 { J ) ,J=1,KQ1) 

READ  3, (TG4( J ) , J=KQ2,KC3) , { TG5 ( K ) , K= 1 , KR 1 ) , ( TG6 { K ) »K=KR2»KR3) 

3, ( ( TC ( L , M ) , L= 1 , LC ) ,M=1»6) 

3, (S(M)  ,M=1,6), (D(M),M  = 1,6),  ( ZL ( M ) ,M=1 ,6)  » (CK(M)  , M= 1 , 6 ) 

3, (CG(M),M=1,6) , (AC(M) ,M=1,6) , (AG(M) ,M=1,6) , (HIM) ,M=1,6) 

3,  ( HD ( M ) » M= 1 , 6 ) , (HG(M) , M= 1 , 6 ) 

3, ( DB (N ) ,N=1,NN) ,( TV (N ) ,N=1,NN) , (DPV(N) ,N=1 ,NN) , (G (N) ,N=1 , NN) 11850401 

1 18504  If 


10 


READ 
READ 
READ 
READ 
6 READ 

DO  55  M=  1 , 6 
I F ( F G ( M ) ) 50 , 5 1 , 50 

51  DO  52  N=1,NN 

52  GSUN(N,M)=0. 

GO  TO  55 

50  DO  54  N=  1 , NN 
q4  CSUN( N, M ) =G ( N ) 

5 
9 
14 


1 185006C 

I 

1 1 8 5 0 0 7 cj 
11851  j 
11851  ! 

1185010 d 
118501  l(j 
1 1 8501 2G 
1185013d 
1 185014C 
1 185015C 
1185  1 
1185017C 
1 185018C 
1 1 8 5 0 1 9( 


CONTINUE 

READ  3ylA,B,C,ZNl,ZN2,ZN37TA,DPA|, 


PB,BS,BL,DGX,DGY,DGZ,DGXl ,DGYllDT 


FORMAT { 4F9 . 4 ) 

READ  14 
GO  TO  300 

13  READ  3, (TGC(K) , K = 1 , KR3 ) 
GO  TO  10 


1185042 C 
1 1850431 
1 1 85044( 
1185045C 
1 185046( 
1 1 85  047( 
1 185048C 
1 185049C 
11851 
11851 
11850521 
1185053C 
1 185054& 


15 

FORMAT ( 43H1 

UNDERGROUND  FALLOUT 

SHELTER 

DATA) 

11850551 

>00 

PRINT  15 

1 1 8 5 0 5 6 ( 

PRINT  16 

1 1 8505  7( 

16 

FORMAT ( 79H0 

ZN 

A B 

C 

V 

G 

1 185058(^ 

1 ZN1  ZN2  ZN3 

DT  ) 

1 1 8 5 059( 

30 

FORMAT ( F9  « 0 , 

F7. 1.2F8 

• 1,F10*0,FI2.0,F7.0, 

3F6 . 2 ) 

PRINT  30 , ZN , 

A , B » C , V » 

G,ZN1,ZN2,ZN3,DT 

1185059; 

PRINT  17 

1 185060( 

17 

FORMAT ( 64H0 

M 

1 2 

3 

4 

5 

118506  1( 

1 6) 

1 1850621 

PRINT  1 8 » ( CK { M ) , M= 1 , 

6 ) 

11850631 

18 

FORMAT ( 10H0 

CK  ( M ) 

6F10.3) 

1185064! 

19 

FORMAT ( 1OH0 

CG  ( M ) 

6F10.3) 

1185065! 

20 

FORMAT ( 10H0 

AC  ( M) 

6F10.3) 

11850661 

21 

FORMATt 10H0 

AG  ( M ) 

6F10-3) 

11850671 

22 

FORMAT ( 10H0 

DIM) 

6F10.3) 

1 1850681 

23 

FORMAT ( 10HO 

Z L ( M ) 

6F10. 3) 

1185069' 

24 

FORMAT ( 10H0 

S ( M ) 

6F10.3) 

1185070! 

5 

FORMAT ( 10H0 

H{  M) 

6F10.3) 

1185071 

<l6 

FORMAT ( 10H0 

HG ( M ) 

6F10.3) 

1185072 

27 

FORMATt 1 OHO 

HD  { M ) 

6F10.3) 

1185073 

PRINT19,  ( CG ( M > ,M=1, 6) 
PRINT20, ( AC(M) , M=l,6) 
PRINT21, (AG(M),M=1,6) 


1185074 

1185075 

1185076 


C°6 


PRINT22 , ( D{ M ) , M=  1 , 6 ) 

1185077! 

PRINT23, ( ZL ( M ) , M= 1 , 6 ) 

1 185078! 

PRINT24, (S(M) , M=l,6) 

1185079! 

PRINT25, ( H{M) , M=l,6) 

1185080! 

PRINT26, { HG ( M ) , M= 1 , 6 ) 

1 18508  1( 

PR  I NT  27 , { HD( M ) , M=l, 6) 

1185082! 

GO  TO  100 

1185083! 

GROUND  TEMPERATURE  INITIALIZATION 

1185084! 

IOO 

IF(WC)  101,102,103 

1185085! 

101 

DO  104  1=1, KP3 

1185086! 

DO  104  J=  1 , KQ3 

1185087! 

DO  104  K=  1 ,KR3 

1185088! 

104 

T G ( I » J , K ) = TGO ( K ) 

1185089! 

GO  TO  129 

1185090! 

102 

DO  105  1=1, KP3 

1185091! 

DO  105  J=  1 , KQ3 

1185092! 

DO  105  K=  1 , KR3 

1185093! 

105 

T G ( I, J,K)  = TGI 

1185094! 

GO  TO  129 

1185095! 

103 

DO  900  1=1, KP3 

1185096! 

DO  900  J = 1 , KQ 3 

1185097! 

DO  900  K=  1 , KR3 

1185098! 

IF(KPl-I)  111,106,106 

1185099! 

106 

IFU-KC1)  111,107,107 

1185100! 

107 

I F ( KQ2- J ) 111,108,108 

1185101! 

108 

IF(K-KRl)  111,109,109 

1185102! 

109 

IF(KR2-K)  111,110,110 

1185103! 

110 

T G ( I,J,K)  =T  G 1 ( I) 

1 185104! 

GO  TO  900 

1185105! 

. .1 

I F ( I — KP2 ) 117,112,112 

1185106! 

112 

IFIJ-KG1)  117,113,113 

1185107! 

113 

I F { KQ2- J ) 117,114,114 

1185108! 

114 

IF(K-KRl)  117,115,115 

1185109! 

115 

I F { KR2— K ) 117,116,116 

1185110! 

116 

TGI  I, J,K)  = TG2  { I) 

1185111! 

GO  TO  900 

1185112! 

117 

IF(KQl-J)  121,118,118 

1185113! 

118 

IF(K-KRl)  121,119,119 

1185114! 

119 

I F ( KR2-K I 121,120,120 

1185115! 

120 

T G ( I,J,K)=TG3(  J) 

1185116! 

GO  TO  900 

1185117! 

121 

IFIJ-KC2)  125,122,122 

1185118! 

122 

IF(K-KRi)  125,123,123 

11851191 

123 

IFIKR2-K)  125,124,124 

1185120! 

124 

TGI  1, J,K)  = TG4l J) 

1185121! 

GO  TO  900 

1185122! 

125 

IF(KRl-K)  127,126,126 

11851231 

126 

TG { I, J,K)=TG5(K) 

1185124! 

GO  TO  900 

1185  125! 

127 

I F { K— KR  2 ) 900,128,  128 

11851261 

128 

T G ( I , J , K ) = TG6 { K ) 

1185127' 

900 

CONTINUE 

1185128' 

1900 

FORM AT l 10F9.5) 

129 

T M=0 • 0 

TQWST=0. 

1185133 

TQWLT=0 • 

1185134 

TTWCT=0 • 

1185135 

T QV  S = 0 • 

1185,136 

T QVL  = 0 • 

1185137 

TQGS=0. 

1185138 

T QGL=0. 

1 1 85  1 39( 

T TCSUN  = 0 « 

1 1 85 140( 

DO  152  M = 1 » 6 

1 185  14 1( 

TQWS ( M ) =0 . 

1185142C 

TQWL(M)=0. 

1185143C 

TTWC ( M ) =0 . 

1185144C 

152 

QSUNT ( M ) =0  . 

1185 145( 

LLC=LC- 1 
ZLC=LLC 

11851 

C END 

OF  INITIALIZATION 

1 1 8 5 146( 

DO  800  N=1,NN 

1185147C 

29 

FORMAT  ( 1 OF  5 . 1 ) 

11851 

301 

DO  28  M= 1 » 6 

28 

PRINT  29,  (TC(L,M),L=1»LC) 

11851 

PRINT  29,  { TGI  I ,KQ,KR  ) , 1 = 1 ,KP1  ) 

11851 

PRINT  29,  (TG(  I , KQ , KR ) , I=KP2,KP3) 

11851 

PRINT  29,  (TG(KP,J,KR) » J = 1 , KQ  1 ) 

11851 

PRINT  29,  (TGIKP, J,KR) ,J=KQ2,KQ3) 

11851 

PRINT  29,  (TGIKP, KQ, K) ,K=l,KRl ) 

11851 

PRINT  29,  (TGIKP, KQ, K) ,K=KR2,KR3) 

11851 

1501 

FORMAT { 16F5 . 1 ) 

PRINT  1501,  TG(  1,  1,1), TG( 1,1, KR3)  ,TG( 1,KQ3,1) , TG ( 1 , KQ3 , KR3 ) 

,TG( KP311851 

1,  1 , 1 ) ,TG(KP3',  1 ,KR3)  , TG  ( KP3  , KQ3 , 1 ) , TG  ( KP3  , KQ3  , KR  3 ) , TG  ( KP  1 , KQ1  , KR  1 ) , 1 1 85  1 

1TGIKP1 ,KQ1,KR2) ,TG(KP1,KQ2,KR1 ) , TG ( KP 1 , KQ2 , KR2 ) , TG( KP2 , KQ1 , 

KR 1 ) , TGI l 85 1 

1 ( KP2 , KQ 1 , KR2 ) , TGI KP2 , KQ2 , KR 1 ) , TG ( KP2 , KQ2 , KR2 ) 

11851 

TM=TM+DTT 

11851 

PSV=GN(DPV(N) ) 

1185148 

WV=0.622*PSV/(P8-PSV) 

1185149 

DTM=0 . 

11851 

1 . j5 

DTM=DTM+DT 
DT  A = 0 . 2 
T A 1 = 40  « 0 

11851 

I R=  1 

1185  1 

405 

QVS=1.08*G*(TA1-TV(N) ) 
QGS=QG(TA1,ZN1,ZN2,ZN3, 1. )+BS 
QWST=0. 

DO  153  M= 1 , 6 

DDI M)=D(M)/ZLC 

QWSIM)=HIM)*SIM)*ITA1— TCI 1,M) ) 

11851 

153 

QWST=QWST+QWS(M) 

GO  TO  I 1153,406) , IR 

11851 

1153 

ZX=QGS-QVS-QWST 
IFIZX)  400,401,402 

11851 

402 

TA1=TA1+DTA 
ZX2=A8SF(ZX) 
GO  TO  405 

400 

T A2=T  A 1-DT  A 
ZX1=AB5F(ZX) 

TA=( TA1*ZX2+TA2*ZX1)/I ZX1+ZX2) 

I R = 2 
TA1=TA 

11851 

GO  TO  405 

11851 

401 

T A=  T A 1 

406 

DDPA=0.2 

DP A 1 = 40  • 0 

I R=  1 

11851 

409 

PSA=GN( DPA1) 

WA=0,622*PSA/( PB-PSA ) 

1185159 

QVL=4780.*G*{ WA-WV) 

118516C 

QWL  T = 0 » 

118516C 

TWCT=0 • 

DO  158  M= 1,6  C"8 

1181 

118| 

157  PS*=GN( TC ( 1.  M ) ) 118516; 

WS=0.622*PSW/ ( PB-PSW ) 118516; 

IF(WA-WS)  154,154,155  1185162 

154  QWL(M)=0.  1185162 

GO  TO  156  118516* 

155  QWL(M)=1061#*HD(M)*S(M)*[WA-WS)  118516* 

1 56  GWLT=GWLT+CWL ( M ) 118516E 

TWC(M)=GWL{M)/1061.  1185165 

158  TWCT=TWCT+TWC(M)  1185165 

QGL  = QG(TA,ZN1,ZN2,ZN3,0.  )+BL  118517C 

GO  TO  { 1158,413) , IR  11851 

1158  ZY=QGL-QVL-QWLT  11851 


IF(ZY)  410,411,412 
412  DPA1=DPA1+CDPA 
Z Y2= ABSF { ZY  ) 

GO  TO  409 

410  DPA2=DPA1— CDPA 
ZY1  = ABSF { ZY  ) 

DP A=  ( DPA1*ZY2+DPA2*ZY1 )/(ZYl  + ZY2) 

I R = 2 11851 


DPA 1=DPA 

GO  TO  409  11851 

411  DPA=DPA1 
413  POA=GN I T A ) 


RHA=1C0.*PSA/P0A 
DO  161  M=  1 , 6 
TGWS(M)=TQWS(M)+QWS(M)*DT 
TQWLIM)=TQWLIM)+QWLIM)*DT 

ZW=CD(M)*(GWS(M)+QWL(M) )/(CK(M)*S(M))+TCl2,M) 

160  TC(1,M)=ZW 

161  TTWCIM)=TTwC(M)+TWC(M)*DT 
TQWST=TQWST+QWST*DT 

I QWLT=TQWLT+QWLT*DT 

TTWCT=TTWCT+TWCT*DT 

TQGS=TGGS+GGS*DT 

TQGL=TGGL+GGL*DT 

TQVS=TGVS^QVS*DT 

TQVL=TGVL+CVL*DT 

C CONCRETE  WALL  TEMPERATURE  DISTRIBUTION 
CO  164  M=  1 , 6 
EIM)=AC(M)*DT/DD(M)**2 
DO  164  L=2 , LLC 

164  TClL,M)  = E{M)*(TC(L-l,M)-*-TClL  + l,M)  + ( 1 . / E ( M ) -2 . ) * TC { L 

I F ( DTM-DTT  ) 1405,1164,1164 

C GROUND  AND  CONCRETE  BOUNDARY  TEMPERATURES 
1164  DO  165  J=KC1,KQ2 
DO  165  K=KR1,KR2 

TGIKP1,J,K)  =FN( CK( 1 ) ,CG( 1 ) , TCI LLC, 1 ) ,TG(KP1-1 , J,K) 

165  T G ( KP2 , J , K ) =FN I CK ( 2 ) , CG ( 2 ) , TC ( LLC  , 2 ) , TG ( KP2+1 , J , K ) 
DO  166  I=KP 1 , KP2 

DO  166  K=KR 1 , KR2 

TGI  I,KG1,K)  =FN(CK( 3), CGI  3) , TCI LLC, 3) ,TG( I ,KQ1-1,K) 

166  TGI  I , KC2 , K ) =FN(CK(4) ,CGI4) ,TC(LLC,4)  ,TG( I ,KQ2+1,K) 
DO  167  I =KP  1 , KP2 

DO  167  J = KQ  1 , KQ2 

TGI  I , J,KR1  ) =FN( CKI 5 ) ,CG( 5) , TC I LLC, 5) , TGI  I , J,KR1~1 ) 

167  TGI  I , J , KR  2 ) =FN I CK ( 6 ) , CG I 6 ) , TC I LLC , 6 ) , TG l I , J , KR2+ 1 ) 
TCILC, 1 )=TGIKP1,KQ,KR) 

TCILC,2)=TGIKP2,KQ,KR) 

TCILC, 3)=TG(KP,KQ1,KR) 

TCILC, 4)=TG(KP,K02t KR) 

TCILC, 5)=TGIKP,KQ,KR1) 

TCILC, 6)=TG(KP,KQ,KR2) 

PRINT  1900, (TCILC, M),M=1, 6)  C_9 


1185175 
1185176 
1185176 
1185177 
1185165 
1185166 
1185177 
1185178 
1185178 
1185179 
1185179 
1185180 
1185  180 
1185181 
1185202 
1 185203 
1185204 
11851 

, M ) ) 1185206 

11851 
1185207 
11851 
1185209 

, DDI 1 ) ,DGX,0. ) 11851 

, DDI2) ,DGX,0. ) 11851 

1185212 

1185213 

, DDI3) ,DGY,0. ) 11851 

, DDI4) ,DGY,0. ) 11851 

1185216 

1185217 

, DDI5) ,DGZ,0. ) 11851 

, DDI6) ,DGZ ,0. ) 11851 
11851 
11851 
11851 
11851 
11851 
11851 


GRQUND  TEMPERATURES  1185220 

KP4  = KP3-1  1185220 


KQ4  = KG3— 

1 

1185220 

KR4-  ? KR3- 

1 

1185220 

00  600  I 

= 2 , K P 4 

1185221 

DO  600  J 

-2 i KQ4 

1185222 

DO  600  K 

= 2 » KR4 

1185223 

IF(KP1 — I ) 

173,  168,  168 

1 185224 

168 

I F ( J-KC1 ) 

173, 169, 169 

1185225 

169. 

IF (KQ2-J  ) 

173,  170, 170 

1185226 

170 

I F ( K— KR  1 ) 

173,  171,  171 

1185227 

171 

IF ( KR2-K  ) 

173,  172,  172 

1185228 

172 

AGG  = AG ( 1 ) 
IF ( I-KP1  ) 

500,600,600 

1185229 

173 

IF{  I-KP2) 

179, 174, 174 

1185231 

174 

I F ( J-KC  1 ) 

179, 175, 175 

1185232 

175 

I F ( KQ2— J ) 

179, 176,  176 

1185233 

176 

IF (K-KR1  ) 

179,177,  177 

1185234 

177 

I F ( KR2— K ) 

179, 178,  178 

1185235 

178 

AGG= AG ( ? ) 
IFIKP2-I  ) 

500,600,600 

1 185236 
1 185238 

179 

IF (KG1-J  ) 

183,  180,  180 

180 

I F ( K — K R 1 ) 

183,  181,  181 

1185239 

181 

I F ( KR2-K  ) 

183,  182,  18? 

1185240 

182 

AGG= AG ( 3 ) 
IF ( J-KG1  ) 

500,590,590 

1185241 

590 

IFl  I-KP1  ) 

500,580,580 

580 

IF  (KP2-I) 

500 , 600,600 

' ° 3 

I F ( J-KG2  ) 

187,  184,  184 

1 185243 

_ _>4 

IF { K— KR  1 ) 

187,  185, 185 

1185244 

185 

I F ( KR 2 — K ) 

187,  186,  186 

1185245 

186 

AGG= AG ( 4 ) 
IF  (KG2-J) 

500,570,570 

1185246 

570 

IF  II-KP1) 

500,560,560 

560 

I F ( K P 2-  I ) 

500,600,600 

187 

IF(KRl-K) 

189, 188,  188 

1185248 

188 

AGG= AG ( 5 ) 
IF  { K— KR  1 ) 

500,550,550 

1 185249 

550 

IF  (I-KP1) 

500, 540, 540 

540 

IF  (J-KQ1) 

500,530,530 

530 

IF  (KQ2-J) 

500,520,520 

520 

I F ( KP2-  I ) 

500,600,600 

189 

1FIK-KR2) 

600,  190,  190 

1 185251 

190 

AGG=  AG  I 6 ) 
I F ( KR2-K  ) 

500, 521, 521 

1 185252 

521 

IF  (I-KP1) 

500,522,5 22 

522 

IF ( J-KG1  ) 

500,523,523 

523 

IF1KG2-J ) 

500,524,524 

524 

IFIKP2-I  ) 

500,600,600 

500 

IF { I-KP1  ) 

502,502,503 

1 185252 

503 

I F ( KP2-  I ) 

502,502,504 

1 185252 

504 

DGX=DGX  1 

1 185252 

502 

I F ( J-KC1) 

501, 501,505 

1 185252 

505 

IFIKQ2-J  ) 

501,501,506 

1 185252 

6 

DGY  = DGY 1 

1 185252 

5T)  1 

£GX  = AGG*DT  r/ ( DGX**2 ) 

11851 

EGY  = AGG*DTT / ( DGY**2  ) 

11851 

EGZ= AGG*DT  T / ( DGZ**2 ) 

11851 

T 1 = TG l 1-1 , 

J , K ) 

1 1 85256 

T 2=  T G ( l + l , 

J , K ) 

1 185257 

C-10 


T 3 = T G ( I»J-1.K)  118  5 

TA=TG(  I , J + l » K ) 1185 

T 5 = T G { I » J » K-  l ) 1185 

T6=TG( I, J,K+1)  1185. 

XX=EGX*(Tl+T2)+EGY*(T3+TA)+EGZ*(T5+r6)  1185. 

YV= 1 .-2 . * ( EGX+EGY+EGZ  ) 1185. 

TGI 1,J,K)=XX+TG( I,J,K)*YY  1185. 

600  CONTINUE  1185, 

C GKGUND  SURFACE  TEMPERATURES  1185. 

DO  191  J = 1 , K Q 3 1185. 

DO  191  K = l , K R 3 1185, 

TG(  1,J,K)=FN(HG( 1 ) , C G I 1) »Db(N)  , 1 G I 2 , J , K ) , 1. ,UGX»QSUNIN, 1 ) ) 1185; 

191  TG(KP3»J,K)=FN( HG { 2 ) , CG I 2 ) , DtJ  I N ) » TG  I K P 3-  1 , J , K ) , 1 • , DG X , Q SUN ( N , 2 ) ) 1185. 

DO  192  1=1, KP3  1185; 

DO  192  K = 1 , KR  3 1185, 

T G ( I , 1,  K ) = FN( HG( 3 ) ,CG( 3 ) , DBI N ) , TGI  I , 2,  K ) , 1 . ,DGY , uSUN I N , 3 ) ) 1185, 

192  TG(I,KG3,K)  = FN(HG(A),CGIA),lJB(N),TG(I,KQ3-l,K),l.  , UG Y , Q SUN  I N , A ) ) 1185, 

DO  193  1=1, KP3  1185; 

DO  193  J=1,KQ3  1185; 

TG(  I , J,  1 ) = FN(  HG(  5 ) ,CG(  5)  ,Dtt(N  ) » TG(  I , J , 2)  , 1 . , DG7. , QSUN  ( N , 5 ) ) 1185; 

193  TGI  I,J,KR3)  = TGI  I,  J,KR3) 


1185; 
1185; 
1185; 
1185; 
1185  j 


AT ( N ) = T A 

1185; 

RH ( N ) =RHA 

1185; 

DP { N ) =CP A 

1185; 

WRITE  OUTPUT 

TAPE  81 ,29,  ( TCI  1, 

M) ,M=l ,6) , TGI  1, 1 , 

1 ) , TGI  1 , 1 ,KR3 ) 

WRITE  OUTPUT 

TAPE  81 , 1900, TM 

195 

FORMAT ( 10H0 

QWSIM)  6F15.0) 

1185; 

196 

FORMAT ( 10H0 

TQWSIM)  6F15.0) 

1185; 

197 

FORMAT { 10H0 

QWL(M)  6F15.0) 

1185; 

198 

FORMAT ( 10H0 

TQWLIM)  6F15.0) 

1185; 

199 

FORMAT { 10H0 

TWC(M)  6F15.0) 

1185; 

200 

FORMAT ( 10H0 

TTWl(M)  6F15.0) 

1185; 

201 

FORMAT ( 10H0 

T C ( 1 , M ) 6F15.1) 

1185; 

202 

FORMAT I 100H0 

1 

2 

3 

1185; 

L A 

5 

6 ) 

lies; 

203 

FORMAT ( 100H1 

N TM 

QVS 

QVL 

QGS1185; 

1 

QGL  GWST 

QWLT  ) 

1185; 

20A 

FORMAT ( 100H 

TWCT 

r qvs 

TQVL 

1185; 

1 TQGS 

T QGL 

TQWST  ) 

11855 

205 

FORMAT ( 1 OOF 

TQWLT 

TTWCT 

TQSUN 

11855 

ITTGSUN 

) 

11855 

206 

FORMAT { 1 I 3, 

1F6.  1.6F15.0) 

1185  5 

207 

FORMAT { 6F  15 

.2  ) 

208 

FORMAT ( 5F  15 

.0  ) 

11855 

302 

PRINT  203 

11855 

PRINT  206,  N 

, TM,QVS,QVL,QGS,QGL ,QWST,QWLT 

11855 

PRINT  20A 

11855 

PRINT  207,  TWCT ,TQVS»TQVL , TOGS 

, TQGL , TQWST 

11855 

PRINT  205 

11855 

PRINT  208,  TQWLT,TTWCT,T-fr5'WrTTQS-W 

11855 

PRINT  202 

11855 

PRINT  195,(QWS(M),M=1,6) 

11855 

PRINT  196, { TQWSIM) ,M=1,6) 

11855 

PRINT  197, IQWLIM) ,M=1,6) 

11852 

PRINT  198, (TQWLIM) ,M=1,6) 

11853 
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PRINT  199, (TWC(M) ,M=1,6) 

1185317 

PRINT  200, (TTWC(M) , M= I , 6 ) 

1185318 

PRINT  201, (TCI  1 , M ) , M= 1 , 6 ) 

1185319 

PRINT  209 

1185320 

FORMAT ( 100H  TV(N)  DPV(N) 

DB  ( N ) 

1185320 

TA  DP A RHA 

1185320 

PRINT  29,  AT(N) ,DP(N) ,RH(N) 

PRINT  2 0 7 , T V ( N ) , DPVIN) , DB ( N ) ,T A, DP A, RHA 

1185320 

DIMENSION  A A ( 2 ) ,EB(2),TN( 170) 

1185340 

A A { 1 ) =0 . 0 

1185340 

AAl 2)=500. 

11851 

8B ( 1 )=100. 

1185340 

BB l 2 ) =50 . 

1185340 

T M = 0 • 

1185340 

DO  90  N= 1 , NN 

1185340 

T M=  TM  + DT  T 

TN ( N ) =TM 

1185340 

CALL  PLOT  ( 3 , NN , TN , A T , NN , TN , DP , 2 , A A , BB ) 

1185340 

CALL  SYSTEM 

1185341 

GOTO  99 

1185341 

ENC 
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Program  M-2  Input  and  Output  Symbols 


1X1 : 

1X2: 

IY1: 

IY2: 

IZ1: 

IZ2: 

IZ3: 

CG: 

AG: 

DT: 

DX: 

DY: 

DZ: 

DTG: 

A,  B, 


A matrix  index  marking  the  boundary  between  the  shelter  air 
space  and  surrounding  earth  (longitudinal  direction) 

A matrix  index  marking  the  end  of  earth  region  along  the 
longitudinal  direction. 

A matrix  index  marking  the  boundary  between  the  shelter  air 
space  and  surrounding  earth  (transverse  direction). 

A matrix  index  marking  the  end  of  earth  region  along  the 
transverse  direction. 

A matrix  index  marking  the  boundary  between  the  shelter  air 
and  roof  earth  region. 

A matrix  index  marking  the  boundary  between  the  shelter  air 
and  floor  earth  region. 

A matrix  index  marking  the  end  of  earth  of  shelter  floor 
region. 

Thermal  conductivity  of  earth,  Btu/hr,ft,  °F. 

2 

Thermal  diffusivity  of  earth,  ft  /hr. 

Finite  difference  time,  hr 

Finite  difference  length  along  the  longitudinal  direction,  ft. 
Finite  difference  length  along  the  transverse  direction,  ft. 
Finite  difference  length  along  the  vertical  axis,  ft. 

Finite  difference  time  for  computing  the  ground  surface  heat 
exchange,  (modification  of  DT),  hr. 

: Internal  dimensions  of  shelter,  ft. 
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CFM:  Shelter  ventilation  air,  cu.ft/min. 

HV  = (60)  (specific  heat  of  air)(density  of  air)  = 1.08  for  standard 
air,  Btu/hr,  °F , CFM. 

WG:  Work  cell  for  ground  surface  temperature  calculation 

if  WG  < 0,  ground  surface  temperature  = outdoor  air  temperature, 
if  WG  > 0,  ground  surface  temperature  is  computed  by  solar  heat, 
convection  heat  and  conduction  heat. 

RUN:  Run  number. 

2 

HW:  Vertical  wall  surface  heat  transfer  coefficients,  Btu/hr,  ft,  °F. 

2 

HR:  Ceiling  surface  heat  transfer  coefficient,  Btu/hr,  ft,  °F. 

2 

HF:  Floor  surface  heat  transfer  coefficient,  Btu/hr,  ft,  °F. 

2 

HG:  Ground  surface  heat  transfer  coefficient,  Btu/hr,  ft,  °F. 

TSW,  TSR,  TSF:  Initial  surface  temperatures  of  wall,  ceiling  and  floor 

respectively,  °F. 

TA:  Shelter  air  temperature,  ®F. 

DPA:  Shelter  air  dew  point  temperature,  °F. 

ZNl:  Number  of  600  Btu/hr  occupants. 

ZN2:  Number  of  400  Btu/hr  occupants. 

Zn3:  Number  of  200  Btu/hr  occupants. 

BS:  Sensible  heat  generated  in  the  shelter  by  non-human  source  .> 

Btu/hr. 

BL:  Latent  heat  generated  in  the  shelter  by  non-human  source, 

Btu/hr . 

PB:  Barometric  pressure,  in.  hg. 

Ll:  if  positive,  initial  ground  temperature  varies  with  depth. 
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L2 : 

If  positive,  outdoor  thermal  environment  is  constant. 

L3: 

If  positive,  ground  temperature  will  be  printed  out  only  at  the 
end  of  total  time  iteration. 

L4: 

If  positive,  ventilation  air  rate  varies  with  time. 

L5: 

If  positive,  number  of  600  Btu/hr  occupants  varies  with  time. 

L6: 

If  positive,  number  of  400  Btu/hr  occupants  varies  with  time. 

L7 : 

If  positive,  number  of  200  Btu/hr  occupants  varies  with  time. 

L8 : 

Strength  of  non-human  heat  source  and  sink  will  vary  with  time. 

NN: 

Total  number  of  time  iterations. 

^PT2 : 

Work  cell  for  two  different  kinds  of  data  input  formats. 

ZLW: 

Lewis  number  for  wall  = 1. 

ZLR: 

Lewis  number  for  roof  = 1. 

ZLF: 

Lewis  number  for  floor  = 1. 

TG*M: 

Constant  initial  earth  temperature,  when  Ll^O,  °F 

TG^(K) : 

Depthwise  variation  of  initial  earth  temperature  when  L1>0,  °F. 

DB^: 

Constant  outdoor  air  dry-bulb  temperature,  °F,  when  L2— 0. 

> 

H 

Constant  ventilation  air  dry-bulb  temperature,  °F,  when  L2^0. 

QSUN^: 

Constant  solar  radiation,  Btu/hr,  ft  , when  L2— 0. 

DPV^: 

Constant  ventilation  air  dew-point  temperature,  °F,  when  L2— 0. 

DB(N) : 

Outdoor  air  dry-bulb  temperature  (time  dependent),  °F. 

TV(N) : 

Ventilation  air  dry-bulb  temperature  (time  dependent),  °F. 

DPV(N) : 

Ventilation  air  dew-point  temperature  (time  dependent),  °F. 

QSUN(N) 

o 

: Solar  irradiation  (time  dependent),  Btu/hr  ft  . 

DN(N) : 

Ventilation  air  rate  (time  dependent),  cfm. 

ZN1N(N) : Number  of  600  Btu/hr  occupants  (time  dependent). 
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ZN2N(N) : Number  of  400  Btu/hr  occupants  (time  dependent). 

QN3N(N) : Number  of  200  Btu/hr  occupants  (time  dependent). 

BSS(N):  Sensible  heat  due  non-human  heat  source  (time  dependent),  Btu/hr. 
BLL(N) : Latent  heat  due  non-human  heat  source  (time  dependent),  Btu/hr. 

9 

QX:  Heat  flux  on  wall  normal  to  longitudinal  axis,  Btu/hr  ft  . 

2 

QY:  Heat  flux  on  wall  normal  to  transverse  axis,  Btu/hr  ft  . 

2 

QYR:  Heat  flux  on  ceiling  surface,  Btu/hr  ft  . 

QZF:  Heat  flux  on  floor  surface,  Btu/hr  ft^. 

TM:  Elapsed  time,  hr. 

STM(N) : Mean  surface  temperature,  °F. 

ST1(N):  Average  temperature  of  the  surface  normal  to  longitudinal  axis,  °F. 

ST2(N):  Average  temperature  of  the  surface  normal  to  transverse  axis,  °F. 

ST3(N):  Average  surface  temperature  of  ceiling,  °F« 

ST4(N):  Average  surface  temperature  of  floor,  °F. 

RHA(N) : Shelter  air  relative  humidity,  %. 

TAA(N) : Shelter  air  dry-bulb  temperature,  °F. 

Y(N) : Elapsed  time  = TM,  hr. 
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C SIMPLIFIED  SYMMETRIC  SHELTER  HEAT  TRANSFER  ANALYSIS  .T.KUSUCA  118510 
C RECTANGULAR  PARALLELEPIPED  SHALLOW  UNDERGROUND  118510 
C S=TCTL  SHELTFP  HEAT  TRANSFER  SURFACE,  SQ.FT  IL35LC 
r CFM=VENTI LATICN  AIR  RATE  ,HV=1.03  SENSIBLE  11851G 


GG 

= TOT  AL 

FEAT  GENERATED  INSIDE 

THE  SHELTER 

113510 

c 

WG 

= WGRK 

CELL,  IF  WG=0. 

.GROUND  SURF ACE =DB ( N ) , IF  NONZERO, CALCULATE 

118510 

c 

IF 

L 1 = 0 

TGC( K )=TGCC 

118510 

c 

IF 

L2  = 0 

DB ( N ) =DBO , T V ( N ) = 

TVO, 

CSU.M  ( N ) = Q S U N 0 

118510 

c 

IF 

L3  = 0 

PRINT  TEMPERATURE 

AT 

EVERY  TIME  INTERVAL 

118510 

c 

IF 

L4=  1 

REAC  CFM  FOR 

E AC  H 

TIME 

PERIOD 

1 1 8 5 1 R 

c 

IF 

L 5 = 1 

REAC  ZN1  FOR 

EACH 

TI*E 

PERIOC 

1 1 8 5 1 C 

c 

IF 

L6=  1 

REAC  ZN2  FOR 

EACH 

TIME 

PERIOD 

118510 

c 

IF 

L7=  1 

REAC  ZN3  FOR 

EACH 

TIME 

PERIOD 

1 1 3 5 I B 

c 

IF  L 8 = 

1 READ  BS  ANC 

BL  FOR  EACH  TIME  PERIOC 

1 1 8 5 1 B 

DIMENSION  STM( 3C0)  ,ST1( 300 ) , ST2 ( 30C ) , ST3 ( 3CC ) , ST4(3C0) 

D I MENS ICN  DPVI300) »Y(30C) ,EB{ 2) ,EC(2) ,RHA( 3GC  ) 

DIMENSION  T ( 20, 20, 4G ) , TGO ( 40) ,DB ( 300 , TV (300) ,QSUN(3C0) 
DIMENSION  TAA( 300 ) , V( 300) ,DN( 3CC ) , ZN1N( 3 CO) , ZN2N( 30C ) , ZN3N ( 30G ) 


DIMENSION  BSS ( 300 ) ,BLL( 300) 

500  FORMAT ( 7 2H  I105LB 

1 ) 1 1 3 5 1 B 

READ  500  1 1 3 5 L 3 

PRINT  500  118513 

2C1  FORMAT ( 1 0p7 . 0 ) 118510 

202  FORM  AT ( 10  I 7 ) 118510 

199  R E AC  202  , 1X1  , 1X2, IY1,IY2,  IZ1, IZ2.IZ3  113510 

READ  201,CG,AG,CT,DX,DY,DZ,DTG,A,6,C  118513 

R E AC  201,  C F M; , H V , WG  , RU  N , HW  , H R » H F , HG  118513 

READ  201,TSW,TSR,TSF,TA,DPA  118513 

READ  201,ZN1,ZN2,ZN3,BS,8L,PB 

RE AC  202 , L 1 , L2 , L3 , NN  » L4 , L5,L6,L7  118510 

R E AC  201,  CPT2 

READ  201,ZLW,ZLR,ZLF 

IF(L1)  203,203,204  118510 

203  READ  201.TGQ0  11851C 

CC  205  K = 1 , I Z 3 118510 

205  TGO ( K ) = TGCG  118510 

GO  TO  206  - 118510 

204  R E AC  2 0 1 *U TGO(K) ,K  = 1,  1Z3)  118510 

206  I F ( L 2 ) 20(7,207,208  1 18510 

207  RE  AC  20 1 , DBG , T VO , QSUNO, CP  VU  113510 

CC  209  N= 1 , NN  116510 

C B ( N ) = DB  0 118510 

T V ( N ) = T V 0 118510 

OPV ( N ) = C PVG  11851 

2 C 9 CSUN  ( N ) = GSUNCI  118510 

GC  TC  210  118510 

208  I F ( C PT 2 ) 401,400,401 


400  RE  AC  402, (CB(N) ,N=1,NN)  , { TV ( N ) , N = 1 , NN ) , (DPV(N) ,N=1 ,NN) , (CSUN(N) ,N  = 


11, NN) 

4C2  FORMAT  (8F9.0) 

GC  TC  210 

4C1  RE AC  201,(CB(N),N=1, NN ) 

REAC  201 , ( TV ( N ) , N = 1 ,NN)  113510 

READ  201, (CPV(N) ,N=1,NN)  11351 

IF  (WG)  61,61,60  1185LB 

60  REAC  201, (GSUN(N) ,N=1,NN)  118518 

61  IF  ( L 4 ) 63,63,62  1 18  5 13 

62  REAC  201 , ( DN( N ) , N=  1 ,NN)  1185  18 

63  IF  (L5)  6 5,65,64  1 1851s’ 
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64  READ  201,(ZN1N(N),N=1,NN) 
*65  IF  (L6J  67 ,67,66 
66  RE AC  20I» (ZN2N[N),N=1,NN) 
" 67“  IF  (TT)  * 6 9,69,68  " 

68  READ  201, (ZN3N(N) ,N=1,NN) 

69  READ  202,  L8 

I F ( L 8 ) 210,210,70 

TO  READ  201, IBSS(N) ,N=1,NN) 
READ  201,  (BLL(N),N=1,NIN) 
210  DU  TO  0~  K=T»  IZ3 
DC  300  J=  1 » I Y 2 
DC  ~ 300  1=1  ,TX2 
I F ( I- I X 1 ) 212,216,216 

I F { I Z 1 — K ) 2 14,216,216 
I FTR-TZZI  *215,716, 21 6' 

GO  TO  300 
TTT » J , K]  =TGO  TK 1 


211 
212 
213 
*214* 

215 

216 
3C0  CONTINUE 

PRINT  * 217, RUN  * 

217  FORM  AT ( 50H 1 SYNMETR ICAL 
PRINT"  218 

218  FCRM AT ( 65H0  S 


1 1 8 5 1 B 

11851B 

1 1 8 5 1 B 

11851b' 

1 1 8 5 1 B 

1 1 8 5 1 B 

1 1 8 5 1 B 

1 1 8 5 1 B 

1 1 8 5 1 B 

118510 

118510 

118510 

118510 

118510 

118510 

1 18  51*0 

118510 

118510 

118510 


SHELTER  HEAT  TRANSFER, RUN 


CFM 


TA 


HV 


F3.0 


HW 


' 1*18  510, 
) 118510 
TT8  510 


1 HG 

SR=A*B 

) 118510 
118513 
riT513 
118513 

SF=SR 

SW=( A+B)*2.*C 

* * ~5=SW+SRTSF  “ 

PRINT  219,S,CFM,TA,HV,HW, HG 

1185*13 

650  FORMAT ( 100H0  WG  HW 

1 ZLR  ZLF  CX 

HR  H F- 

DY  DZ 

ZJW 

) 

PRINT  650 

651  FORMAT! 10F10.3) 

PRINT  651fWGfHW|HRfHFfZLWfZLRfZLFfDXfDY  tDZ 
219  FORM  AT { 6 F 10 • 3 ) 

118510’ 

PRINT 220 

118510 

220  FORMAT ( 55H0  DB ( N ) TV(N) 

QSUN(N) 

1118510 

DC  222  N=1 ,NN 

118510 

222  PRINT  221,  DB { N ) ,TV ( N) , QSUN ( N ) 

118510 

22T'*in:RMTTT*TFT0T2T' 

H85T0 

READ  202 , I I , J J, KK 

C END  OF  DATA  INPU1 

* 118510 

KX= I X2~l 

118510 

**KXT=TXT- 1 * 

1T8510 

KY=I Y2-1 

118510 

KYY=1Y1-1 

' * * * I185TO 

KZ=IZ3-1 

118510 

KZ 1= 1 Z 1~ 1 

TT  851*0 

K Z 2= I Z 2- 1 

118510 

tX=AG*UT / ( UX*UX ) 

118  5T0 

EY=AG*DT/(DY*DY) 

118510 

EZ=AG*DT/{  t)L*DZ) 

***** IT8T5TO 

H1=HW*DX/CG 

118513 

H2=  PW*UYTC*G 

118513 

H3  =HR*CZ/CG 

118513 

H4  = H F * DZ / C G 

118513 

HGU=HG*CZ/CG 

11.8510 

R0=  1 • ~-2  • * ( EX+EY+fcZ  ) 

TT55T0 

RX=RC-2.*EX»H1 

K Y = K U— 2 • * t Y *nc 
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RZR=R0-2.*EZ*H3 

RZF=R0-2.*EZ*K4  

RG=l.-2.*( EZ)*( l.+HGU)*CTG/DT  118510 


( TM  = 0.  ” *■"* 118510 

W'  END  OF  INITIALIZATION  118510 

PRINT  223  " ' 118510 

22 3 FORMAT ( 100H0  RO  RX  RY  RZR  RZF  118510 

1 HI  HZ  H3  H4  ) 

PRINT  229,R0,RX,RY,RZR,RZF,H1,H2,H3,H4 

329  PRINT  228  - - - ■ - * 

228  FORMAT ( 105H  CG  AG  HGU  DX  DY  118510 

1 DZ  CT  EX  EY  EZ  ) 118510 


PRINT  229, CG, AG, HGU  , DX , DY , DZ , CT , EX , E Y , EZ 

229  FORMAT  ( 10FI0 . 31 118510 

234  FORMAT ( 20F6 . 1 ) 118510 

YY= I Y1 

Z Z=  I Z2- 1 Z 1 ' ~ ‘ ’ 

AA=CX*XX*2. 


BE=C Y*YY*2  . 

CC=CZ*ZZ 

SX  = A A*CC  — 

SY  = B B*CC 

SZR=AA*BB 

SZF=SZR 

612  FORMAT C 50  HI  QX  " GY  QZR 

PRINT  612 

wa=w(dpat~ 


WW=k(TSW) 

• WR  = MTSR) 

WF  = W ( TSF  ) 

C START  OF  TIME  ITERATION 
DO  907  N= 1 , NN 

IF”  ( LA  ) 51,51, 5C — 

50  CFM  = DN(N) 

51  IF  CL5)  53,53,52 
52  ZN1  = ZN1N(N) 

53  IF  { L6 ) ”55,55,54  ” 

54  ZN2  = ZN2N ( N ) 

55  IF  (17  T 73,71, 56  

56  ZN3  = ZN3N ( N ) 

71  IFIL8)  57,57,72 
72  BS=BSS ( N ) 

HI  = B L L ( NT 


57  TM=TM+DT 


TTT  VTA- WW) 501  ,*5771 , 602" 

602  ZZW=106i.*( WA-WW)/{ TA-TSW+O . OOOC 1 ) /O. 243 
7ZW=ZZW/ZLW 


GO  TO  603 

6 Cl  ZZW=0~.  

603  IF(VlA-WR  1604,604,605 
“ 6 05  7ZR=T0  61.  * ( WA-WR  77TT  A-TSR  TUZOOCOIT  Z07Z43 

ZZR=ZZR/ZLR 

GO  TO  606""  1 

•604  Z ZR=0 . 

606  IF  ( WA-WF  ) 607 , 607 , 608 " 

608  ZZF=1061.*(WA-WF)/(T  A-T  SF+0.Q0CC1 1/0.243 

ZZF=ZZF/3LF 

GO  TO  609 

607  ZZF=0 . 


QZF 


TM  ) 118513 

118513 

rr85T3- 

118513 

118513 

118513 

1185l0 

118510 

178511 

11851B 
IT851B 
1 1 85 1 B 
IT85I8 
11851B 

1T85TB 

1 1 8 5 1 B 
rr85nr 
118518 
TT85T8 
.118510 

rT8"5T3 

118513 


118513 


118513 

118513' 

"TT85T3 

113513 

118513 


118513 
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60  TSM=(TSW*Sto+TSR*SR+TSF*SF)/(SW+SR+SF)  118513 

PW=  T A 
PR  = T A 

PF-TA  • 

H X = H 1 * ( l.+ZZW)  118513 

H Y = F 2 * ( l.  + ZZW)  " " 118513 

HZR=H3*( l.+ZZR)  118513 

HZF=H4*(  l.+ZZF ) 118513 

QUU=CSUN(N)*DZ/CG  118510 


QS=GG(TA,ZN1 »ZN2»ZN3, l.T+BS 
TA1=TA+1. 

CS1=GG(TA1,2N1,ZNZ ,21X3,  1 .1  + B S 

Y1=CS1-GS 

Y2  = CS*TAT-GSr*TA 

CL=CG(TA,ZN1,ZN2,ZN3,0. )+BL 

WV=WIDPV{NM 

RX=R0-2.»EX*HX 

RY=RD-2 . *EY*HY 

RZR=RO-2.*EZ*HZR 


RZF=R0-2 • *EZ*HZF 
702  IF(L3)231»231,845 

84  5 IF(N-NN)  245,231,231 

231  PRINT  232 

232  FORM AT { 5 0 H 1 EARTH  TEMPERATURE  ON 
DO  235  K= 1 , I Z 3 

235  PRINT  23  V,  ITT  I',  I ,KT  , T=  1 , 1X2 ) 
PRINT  236 

236  FORMAT! 50HX"  EARTH  TEMPERATURE  ON 
DO  237  K=  1 , I Z 3 

z-x  PLANE  AT  J=1 
2-X-  PLANE  at  J=  j j 

118511 
) 1185 11 
118511 
TI8511 
118511 
~ 1TT85X1 
118511 

237  PRINT  2 34  » I TIT,  J X »K '}  » I = 1 » I X 27 
PRINT  238 

238  FORMAT! 5 OH 1 EARTH  TEMPERATURE  ON 

DO  239  K=1,IZ3 

■“  239  PR  TNT  234,TT(  1,  J,KX,  J=1  ,T  Y2T 
PRINT  240 

"240“  F0RHATC5WX  T:  ARTTTT  E M PE  RATURE""ETN 
DO  241  K=  1 , I Z 3 

Z-Y  PLANE  AT  1=1 
Z-Y  PLANE  AT  1=11 

118511 
118511 
"IT1 85 11 
118511 

'118511 

118511 

1118511 

118511 

241  PRINT  234,  (Till, J,K),  J — 1 , IY2) 
PRINT  242 

"242  FORMAT (5 OKI  "EARTH  TEMPERATURE  ON 
DO  243  J=  1 , I Y 2 

' X-Y" PLANE AT'  “K=“KK 

118511 

118511 

11X3511 

118511 

243  PRINT  234, TTII, J»KK) ,1=1, IX2T 
C EARTH  TEMPERATURE  BOUNDARIES 

118511 

118511 

245  DO  103  1=1, 1X2 

DO  103  J= 1 , I Y2 

118511 

118511 

IFIWGl  101,101,102  ’ 

1C1  T(  I, J, 1 )=DB(N) 

118511 

118511 

GO  TO  1C3 
102  TG  = 0. 

110511 

“802  T ( I7J7T)  =2  . *EZ* (TIT,  J, 2H-HGU*DB { N } + QUU ) *DTG7T)T  +TTI7T,  XX*RG 
il1  (TG-CT)  803,103,103 

“"803  TG=TG+DTG~  ' ~ 


GO  TO  802 

~nn~TTT7XriTJT=TGTT  T IZ3)  118511 

TGO  ( 1 )=T  ( 1X2, IY2, 1) 

UNDISTURBED  EARTH  TEMP ERATURE 

DO  104  K=  2 » K Z 118511 

IU4~TU0TR ) = E2*  ( T GCTK'“XX+  TGOTK+ 1 TT+TG  0 ( 10^117=27*^77 1X8511 

DO  105  K=2,KZ 

DO  105  J=  1 , IY2  118511 
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1C5  T (1X2,  J,K)=TGO(K)  118511 

DC  r06  K=2,K2  ' ' 

DC  106  1=1, IY2 

TC6  T'(  I , TY2 , K I=TG0 ( K ) " — 


DC  900  K=2 , KZ 

DC  900  J=1,KY 

DC  900  1=1, KX 

P2=T1  T>  r,^j,  K ) 

P4  = T ( I » J+l , K ) 

P5  = T f I , J , K-l ) 

P6=T { I , J » K + 1 ) 

I F { I-  1 ) 1,1,2 
1 P 1=P2 

G C TD  3 


2 P 1=T { I - 1 , J,K) 
“3"  IFTJ-TT  4,4,5 
4 P3=P4 


- GC  T Q g 

5 P3=T ( I , J- 1 , K ) 

6 P0=T (T,J , K ) 

C*  END  CF  PREPARATION 


TFrr^TxiT-  7 , ii  9 i i ~ 


7 I F ( J-IY1 ) 8,11,11 

8 IF ( IZ1-K  ) 9, 11,11 

9 I F ( K- I Z2 ) 10, 11,  11 

10  GC  TO  900 

11  I F ( K-I  Z 1 ) 8 00  » 12, 19 
"T2TTFTT-I X O 13 , 15 , 800' 


13  IF(J-IYl)  14,18,800 

14  P0=EX* ( P1  + P2 ) + EYMP3  + P4 ) + 2. *EZ* ( P5+HZR+PR ) +P 0 * R Z R 
GC  TO  900 


1 51 F { J- I Y 1 ) 16,17, 800 

16  P0=( EX/3. )* ( 2.*P1+4.*P2 ) + EY* { P3+P4 ) + ( EZ/3 . )*(4.*P5+2.*P6)+P0*R0 

to  to'  900"  • ~ 


118511 
118511 
118511 
118511 
118511 
118511 
118511 
118511 
118511 
1 18  5 1 1 
118511 
118511 
118511 
118511 
118511 
118511 
11851L 
118511 
11851  1 
118511 
118511 
118511 
118511 
118511 
118511 
118513 
118511 
118511 
118511 
118511 


17  P0=  ( EX  / 7 • )*( 6.*P1+8.*P2 )+ (EY/7.  )*(6.*P3+8.*P4)+(EZ/7.  )*(8.*P5+6.*P1185ll 

16 ) +P0*RC  ‘ ' “ " 118511 

GO  TO  900  118511 

18  PO=EX *( Pl  + P 2 ) + ( E Y/3 . ) * ( 2 . *P  3+4. *P4 )+{ EZ/3. ) * ( 4 . * P 5 + 2 . * P 6 ) + PO * R C 118  511 

GO  TO  900  118511 

T9~TfT1Z2-KT  800 , 26 t 29~ 

20  IF(I-IXl)  24,21,  800  118511 

2TTFTJ-T  YTJ  22,2 3,  ~ SCO  " 1 185 1 1 


22  P0=EY*(P3+P4)+EZ*(P5+P6)+2.*EX*(P2+HX  *PW)+PO*RX  118513 


23 

TO  TO  900 

P0=( EX/3. )*(2.*P1+4.*P2 )+ (EY/3. ) * ( 2 . *P3+4 . *P4 ) + EZ * ( P 5 + P6 ) +PO*RG 

118511 

118511 

24 

GTTTfTTOcr 

I F ( J-IY1 )900,25,  800 

118511 

118511 

25 

~PT5=EX*(  P1  + P2)+EZ>(P5  + P6  ) + 2.*EY*{  P4+HY  »PWT+PO*RY 
GO  TO  900 

•118513 

118511 

26 

27 

I F { I - 1 X D 27  ,T0 , 800 
IF(J-IYl)  28,29,800 

118511 

118511 

28  PO=EX*( P1+P2)+EY*(P3+P4)+2.*EZ*( P6+HZF*PF )+PO*RZF 
GO  TO  900 

118513 

118511 

29 

PO=EX*( Pl+P2)+(EY/3. )*( 2.*P3+4.*P4)+( EZ/3. ) * ( 2 . *P 5+4. *P6 } + PO*RO 
GC  TO  900 

118511 

118511 

9 

30 

31 

IF(J-IYl)  31,32,800 

P0=( EX/ 3. )*{ 2.*P1+4.*P2 )+EY*( P3+P4 )+( EZ/3. )*( 2 . *P 5+4. *P6 ) + PO*RO 

118511 

118511 

32 

GO  TO  900 

P0=( EX/7. )*{ 6.*P1+8.*P2 )+(EY/7. )*( 6 .*P3+8 . *P4 ) + ( E Z/7. )*(8.*P6+6. 

118511 

•P113511 

1 5 ) +PO*RO 

118511 
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118511 

118511 


GC  TO  900 

800  PC=EX*tPl+P2)+EY*(P3+P4)+EZ*( P5+P6)+PQ*R0 
9C0  T { I , J , K ) =P0 

910  ZT=0 . 

ZN  = 0 . 

DC  911  K=IZ1»IZ2 

DC  911  J=  1 , I Y 1 

ZT=ZT+T( 1X1, J,K) 

911  ZN=ZN+ 1 . 

TSX=  ZT/ZN 
Z T = 0 . 

ZN  = 0 . 

DC  912  K=IZ1,IZ2 

DC  912  T=T7TX1 

Z T=  ZT  + 1 ( I , I Y 1 1 K ) 

912  ZN=7N+T“.  

TSY=ZT/ZN 

ZK=0  - 

ZT  = 0. 

DO  913  1 = 1. 1 X 1 

DO  913  J=1,IY1 

7T=zT+rn;j, IZ1)  - ~ * 

913  ZN=ZN+1. 

TSZR=ZT/ZN 
ZN  = 0 . 

ZT=0  . 

DO  914  1 = 1, 1X1 

TKJ—  914  JT=  1 , 1 Y 1 ~ 

ZT=ZT+T ( I , J, IZ2) 

~9“14P  ZN=ZN  + 1.  “ 

TSZF=ZT/ZN 


11851 
118  51 
"11851 
11851 
11851 
11851 
11851 
11851 
11851 
11851 
118  51 
11851 
11 8 5 1 
11851 
11851 
11851 
' 118 51 
11851 
11851 
11851 
11851 
11851 

1185  r 

11851 
1 1851 
11851 
1 18  51 
11851 


GX=HW*( l.+ZZW) * ( PW-TSX1 
QY=HW*( 1 .+ZZW )*{ PW-TSY) 
QZ'K^Tf  FT*TT . + Z Z RT*TP'R-TSZR  T “ 
QZF=HF*U.+ZZF)*f  PF-TSZF) 
TSW=  (T3X  *SX+T  SYl^TTrSXTTYT 
TSF=TSZF 

TSR=TSZR  ‘ _ ~ 


118513* 
118513* 
118513* 
118  513*- 
1 1 8 513* 
1 185 .13* 
118513* 


TSM=( TSW*SW+TSR*SR+TSF*SF  )/ ( SW+SR  + SF)  118513' 

TA=T1 . 0 8 *C  F M*  T V ( N ) + Y2+TIW*  TS  VI*  SV+HR*T$R*SR  +T)F*tSF*  SFf  7T1 . 08*CFM-Y  1 + 1 18  5 1 3 f 


lHVi*SW  + HR*SR  + HF*SF  ) 118513' 

~ G L = C G I T7,  m 1 , 7N 2 »~ZN 3,Tr.T+TJL  1185I31 

WW= W { TSW ) 118513' 

ftR  = WTT5 R J II  85TT‘ 

WF  = to ( TSF ) 118513? 

"W=prvf/zi  w 1 — 

HR=HR/ ZLR 

HF=HF/ZLF  


WA=(4780.*CFM*WV+GL+( 1061  . ) * ( WW* SW *HW+WR*SR*HR +WF* SF*HF ) /O . 243 ) / ( 41 1 851 3! 


T780 . *CFM+  ( 1061 . ) * f SW*TW+SR * HR+5T  * H F 1 70 . ^43 T ' 1 1 8 5IT? 

HW=HW*ZL W 

KR'=FTR*ZL  R 

HF=HF*ZLF 

WS=MTSM)  IT55T3T 

IF(WA-WS)  610,610,611  118513( 

6 1 0 ' W A = 1 4 7 8 0 . * C F M*WV+  Q L J 7r47TTU7*FCTMl  lT85lT< 

611  PSA  = WA*PB/( WA+0.622)  . 118513< 

r h~atn  t =raavi  PTATG^rm 

513  FORMAT ( 5F10.2 ) 118513* 

P RT  NT  6 11  ,QX  7QY7  CZRTGZF,"Tft IT  F5T3^ . 
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STM(N)=TSM 
ST1 (N)=TSX 
ST2(N)=TSY 
ST3(N)=TSZR 
ST4(N)=TSZF 
V { N ) = T M 
Y ( N ) =TM 
9C7  T A A ( N ) = T A 
PRINT  950 

950  FORMAT ( I00H1  TM  STM  ST1  ST2  ST3 

1 S T A RHA  TA  A ) 

DC  951  N= 1 » NN 

951  PRINT  954,Y(  N ) ,STM(N  ),  STKN)  ,ST2(N)_,  ST3 ( N ) , S T4 ( N ) , RH A ( N ) , T AA ( N ) 

954  FORMAT ( 8F 10 . 2 ) 

E B ( 1 ) =0 • ' 

EB ( 2 ) = 500 . 

EC ( 1 )=1C0. 

EC ( 2 ) = 40  . 

CALL  PLOT ( 3 » NN*  Y » STM ,NN  » V , T A A , 2 * EB » EC ) 

9C1  CONTINUE  “ 118512' 

CALL  SYSTEM  118512' 

GO  TO  199  ‘ ""  - - ---  - 118512' 

END 
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Program  M-4  Input  and  Output  Symbols 


NN: 

Total  number  of  time  iterations. 

II: 

Number  of  ventilation  air  inlets  where  temperatures  were  observed 

JJ: 

Number  of  shelter  air  temperature  stations. 

KK: 

Total  number  of  shelter  inner  surface  temperature  stations. 

KWC : 

If  —0,  Ft.  Belvoir  1000-man  shelter  analysis. 
If  >0,  Ft.  Belvoir  200-man  shelter  analysis. 

WC: 

If  <0,  constant  earth  temperature,  TG$,  is  read  in. 

If  >0,  six  directional  earth  temperature  profiles,  TG(I,M), 
will  be  read  in. 

AI: 

If  >0,  air  conditioning  calculations  are  included. 

WC1: 

If  —0,  cfm  varies  with  time. 

WC2 : 

If  <0,  number  of  occupants  varies  with  time. 

TG$$: 

Constant  initial  earth  temperature,  °F . 

TG  ( I , M) 

:Six  directional  initial  earth  temperature  profiles,  °F. 

TC^(M)  : 

Initial  inner  wall  temperatures,  °F. 

TC(L,M) 

: Initial  temperature  profile  within  the  shelter  inner  walls,  °F. 

S(K): 

Shelter  inner  surface  area,  ft  . 

DD(K) : 

Finite  difference  length  for  inner  wall,  ft. 

DG(K) : 

Finite  difference  length  for  earth  region,  ft. 

D(K)  : 

Thickness  of  the  inner  wall,  ft. 

ZL(K) : 

Thickness  of  the  earth  block,  ft. 

AG(K)  : 

Thermal  diffusivity  of  earth,  ft^/hr. 

AC  (K)  : 

Thermal  diffusivity  of  inner  wall,  ft  /hr. 
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CK(K)  : 
CG(K) : 
H(K)  : 
HD(K)  : 
DT: 

DTT: 

GS: 

GL: 

GA; 

FR: 

CF: 

DATT: 

TW1: 

XDB(N) : 
TV(N) : 
DPV(N) : 
EG(K) : 

E(K): 

EC(K)  : 
TM: 


Thermal  conductivity  of  inner  wall,  Btu/hr  ft,  °F. 

Thermal  conductivity  of  earth  block,  Btu/hr,  ft,  °F. 

Heat  transfer  coefficient  of  the  inner  surface,  Btu/hr,  ft  , °F. 

o 

Vapor  transfer  coefficient  of  the  inner  surface,  lb/hr,  ft  (lb/lb) . 
Finite  difference  time  for  inner  wall  temperature  calculation,  hr. 
Finite  difference  times  for  earth  temperature  calculation,  hr. 

Per  capita  sensible  heat  generation  from  non-human  heat  source, 
Btu/hr,  person. 

Per  capita  latent  heat  generation  from  non-human  heat  source, 

Btu/hr,  person. 

Recirculation  air  rate  for  conditioning,  cu  ft/min. 

Coil  effectiveness  of  the  air  conditioner  (dimensionless) 0 
Coil  contact  factor  of  the  air  conditioner  (dimensionless). 

Coil  leaving  air  temperature  rise  due  to  fan  heat,  °F. 

Entering  temperature  of  the  coolant  to  the  air  conditioner 
coil,  °F . 

Outdoor  air  dry-bulb  temperature,  °F. 

Ventilation  air  dry -bulb  temperature,  °F. 

Ventilation  air  dew-point  temperature,  °F. 

Finite  difference  stability  modulus  for  earth  temperature 
calculation,  which  should  not  exceed  0.5. 

Finite  difference  stability  modulus  for  shelter  inner-wall 
temperature  calculation,  which  should  not  exceed  0.5. 

E(K)0 

Elapsed  time  (hr) . 
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QVS : 

Sensible  heat  carried  out  by  ventilation  air,  Btu/hr0 

QVL: 

Latent  heat  carried  out  by  ventilation  air,  Btu/hr0 

QGS: 

Sensible  heat  generated  in  the  shelter,  Btu/hr. 

QGL: 

Latent  heat  generated  in  the  shelter. 

QWST: 

Sensible  heat  absorbed  by  shelter  surface,  Btu/hr. 

QWLT: 

Latent  heat  absorbed  by  shelter  surface,  Btu/hr. 

TWCT: 

Condensate  collected  in  the  shelter,  lb/hr. 

QAS: 

Sensible  heat  absorbed  by  the  air  conditioner,  Btu/hr. 

QAL: 

Latent  heat  absorbed  by  the  air  conditioner,  Btu/hr. 

TA: 

Shelter  air  dry-bulb  temperature,  °F. 

DP  A: 

Shelter  air  dew-point  temperature,  °F. 

WBA: 

Shelter  air  wet-bulb  temperature,  °F. 

RHA: 

Shelter  air  relative  humidity,  %. 

EFS1 : 

Shelter  air  effective  temperature,  °F. 

TSM:  Average  shelter  inner  surface  temperature,  °F. 


QWS  (K)  : 

Sensible  heat  transferred  to  the  Kth  surface,  Btu/hr. 

QWL(K)  : 

Latent  heat  transferred  to  the  Kth  surface,  Btu/hr. 

TWC(K)  : 

Condensate  collected  onto  the  Kth  surface,  Ib/hr. 

FLX(K)  : 

Heat  flux  of  the  Kth  surface,  Btu/hr,  ft  . 

TS(K) : 

Average  temperature  of  the  Kth  surface,  °F. 

TC (L,K) 

: Temperature  distribution  in  the  Kth  inner  wall,  °F. 

TG( I, K) : Temperature  distribution  in  the  Kth  block  of  the  earth. 

TIME(N) : Elapsed  time  for  the  observed  date,  hr. 

CFM(N) : Ventilation  air  rate  used  in  the  experiment  (time  dependent) , cfm 
ZN(N):  Number  of  occupants  used  in  the  experiments  (time  dependent). 


026 


XI:  = DBV(N)  - Observed  ventilation  air  dry-bulb  temperature,  °F . 

X2:  = WBV(N)  - Observed  ventilation  air  wet-bulb  temperature,  °F. 

DBS:  Observed  shelter  dry-bulb  temperature. 

WBS:  Observed  shelter  wet-bulb  temperature. 

ETS1:  Shelter  effective  temperature  based  upon  DBS  and  WBS. 

Subscripts : 

MwK  = Type  of  wall  exposures 

1.  North  wall 

2.  South  wall 

3.  East  wall 

4.  West  wall 

5.  Floor 

6.  Roof 
N = Time 

L = Spatial  matrix  index  for  inner  wall 

I = Spatial  matrix  index  for  earth  region 

Subroutines : 

(1)  PREP:  This  subroutine  is  designed  to  accept  observed  data  of 

prototype  shelter  and  yield  the  psychrometrical ly  processed  mean 
shelter  conditions  in  terms  of  temperature,  humidity  and  heat  flux. 
Following  are  the  input  data  specific  to  this  subroutine. 

KWS:  The  first  temperature  entry  for  XWS  table. 

NWS:  The  total  number  of  temperature  entries  for  XWS  table. 

XWS(I):  Saturated  air  humidity  ratio  table  of  Goff  and  Gratch,  lb/lb. 

KHA:  The  first  temperature  entry  for  XHA  table. 
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NBA: 
XHA(I)  : 
KHS : 

NHS : 
XHS(I)  : 
KFS: 
NFS: 
XFS(I)  : 

KHW: 

NHW: 

XHW(I)  : 

KVA: 

NVA: 

XVA(I)  : 

KVS: 

NVS : 
XVS(I) : 
KQS: 
NQS: 
XQS(I) : 
EDB(I): 


EWB(I) : 


The  total  number  of  temperature  entries  for  XHA  table. 

Dry  air  enthalpy  table  of  Goff  and  Gratch,  Btu/lb. 

The  first  temperature  entry  for  moisture  XHS  table,  °F. 

The  total  number  of  temperature  entries  for  XHS  table. 

Moisture  saturated  air  enthalpy  table  of  Goff  and  Gratch,  Btu/lb 
The  first  temperature  entry  for  the  XFS  table. 

The  total  number  of  temperature  entries  for  XFS  table. 

Factors  relating  the  degree  of  saturation  to  the  relative 
humidity. 

The  first  temperature  entry  for  the  XHW  table. 

J 

The  total  number  of  temperature  entries  for  XHW  table. 

Liquid  water  enthalpy  table  of  Goff  and  Gratch,  Btu/lb. 

The  first  temperature  entry  for  the  XVA  table. 

The  total  number  of  temperature  entries  for  XVA  table. 

Dry  air  volume  table  of  Goff  and  Gratch,  cu  ft/lb. 

The  first  temperature  entry  for  the  XVS  table. 

The  total  number  of  temperature  entries  for  XVS  table. 

Moisture  saturated  air  volume  table  of  Goff  and  Gratch,  cu  ft/lb 
The  first  temperature  entry  for  the  XQS  table. 

The  total  number  of  temperature  entries  for  XQS  table. 

Per  capita  sensible  heat  table  of  human  body,  Btu/hr. 

Dry-bulb  temperature  entries  to  ASHRAE  effective  temperature 
table,  °F . 

Wet-bulb  temperature  entries  to  ASHRAE  effective  temperature 
table,  °F . 
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ET(I, J) jASHRAE  effective  temperature  table. 

(2)  STRCON : Subroutine  to  compute  adiabatic  shelter  air  conditions  by 

accepting  ventilation  air  dry-  and  wet-bulb  temperatures  and  ventilation 
air  flow  rate. 

(3)  AIRC^N:  Subroutine  to  compute  sensible  and  latent  heat  absorption 

by  accepting  the  following  input  data. 

GA:  Recirculation  air  rate,  cfm. 

CFM:  Ventilation  air  rate,  cfm. 

DBS:  Shelter  air  dry-bulb  temperature,  °F . 

WBS:  Shelter  air  wet-bulb  temperature,  °F. 

DBV:  Ventilation  air  dry-bulb  temperature,  °F. 

WBV:  Ventilation  air  wet-bulb  temperature,  °F. 

TW1:  Inlet  temperature  of  the  air  conditioner  coolant,  °F. 

E:  Effectiveness  of  the  air  cooling  coil. 

CF:  Contact  factor  of  the  air  cooling  coil. 

DTA:  Dry-bulb  temperature  rise  of  the  air  conditioner  outlet  air 

due  to  the  fan  heat,  °F. 

(4)  FN:  Boundary  temperature  calculation  subroutine  by  a finite 

difference  formula. 

(5)  EFT I:  Subroutine  for  computing  effective  temperature  by  using  the 

input  data  of  dry-  and  wet-bulb  temperatures. 

(6)  DBWBH:  Subroutine  for  computing  enthalpy  of  moist  air  using  dry- 

and  wet-bulb  temperatures. 

(7)  DBWBRH:  Subroutine  for  calculating  relative  humidity  for  given 

dry-  and  wet-bulb  temperatures  and  barometric  pressure,  PB  (in.  Hg) , 
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(8)  DBWWBH:  Subroutine  to  calculate  wet-bulb  temperatures  and  enthalpy 


of  moist  air  based  upon  given  values  of  the  dry-bulb  temperature  and 
humidity  ratio. 

(9)  PV:  Subroutine  for  calculating  the  vapor  pressure  in  inches  of  Hg 

of  saturated  moist  air  at  given  temperatures. 

(10)  DBWBDP : Subroutine  for  calculating  the  dew-point  temperature  and 

humidity  ratio  of  moist  air  for  given  dry-  and  wet-bulb  temperatures. 

(11)  DBWPWB:  Subroutine  for  calculating  the  wet-bulb  temperature  and 

humidity  ratio  when  dry-bulb  temperatures  and  dew-point  temperatures 
are  given. 

(12)  TBLU : Table  look-up  and  linear  interpolation  subroutine  when 

the  independent  variable  is  incremented  by  unity. 

(13)  XTBLU : Table  look-up  and  linear  interpretation  subroutine  when 

the  incrementation  of  the  independent  variable  is  non-unity. 
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C ONE  DIMENSIONAL  SHELTER  II851003 

C SUBROUTINE  TO  BE  USED  IN  THE  PROGRAM  11851004 

C PREP(DT , II , JJ,KK,KWC,NN)  11851M03 

C STRCON (DB,WB,G)  ADIABATIC  SHELTER  AIR  CONDITION  CALCULATION 

C AIRCCN(GA,CFM,DBS,WBS,DBV,WBV»TWI ,E,CF,  QAS,QAL,DTA)  11851A02 

C FUNCTION  FN 

C EFT1 (DB ,W8 »EFX) 


c 

DBWBH ( DB , WB , H ) 

GIVEN 

DB 

AND  WB  , COMPUTE 

H 

11851012 

c 

DBWBRH(DB,WB,PB,RH) 

GIVEN 

DB, 

WB , PB  , COMPUTE 

RH 

11851014 

c 

DBWWBHC  DB,W,WB,H) 

GIVEN 

DB 

AND  W , COMPUTE 

WB 

AND 

H 

11851011 

c 

FUNCTION  PV ( X , PB ) 

CALCULATE 

VAPOR  PRESSURE 

FOR 

GIVEN  DB 

AND  PB 11851018 

c 

DBWBDPCDB,WB,DP,W) 

GIVEN 

DB 

AND  WB  , COMPUTE 

DP 

AND 

W 

11851010 

c 

DBDPWB ( DB»DP,WB,  W) 

GIVEN 

DB 

AND  DP,  COMPUTE 

WB 

11851013 

c 

TBLU  C X, KX, V , Y) 

KX=  FIRST  V 

ARI ABLE 

11851015 

XTBLLMX,  DX,KX,LX,V,Y)  , 


TEMP.  DATA 
TEMP.  DATA 

MAN  INPUT 


11851016 

11851017 

11851005 

11851006 

11851007 

11851008 

11851009 

11851008 

11851M12 

11851004 

11851005 

11851006 

11851007 


V=TABLE 
Y=V ALUE  FOR  X 
DX= VAR  I ABLE  INCREMENT 
KX=FI RST  VARIABLE 
LX=NUMBEROF  ENTRIES 

V =T  ABLE 

Y = VALU E OBTAINED  FOR  X 

1=N0RTH » 2=  SOUTH*  3= EAST  ,4= WE ST »5=FL00R»6=RQ0F»J  = R00F  MATRIX 
AI =0  - NO  AIR  CONDI  TIDING*  +=AIR  CONDITIONING 
WC=0-»  CONSTANT  T GO  ,+=PROFlLE 

II  =NUMBER  OF  VENTILATION  TEMP.  DATA 
JJ  =NUMBER  OF  SHELTER 
KK  =NUMBER  OF  WALL 
KWC  =*  200MAN  INPUT  ,0-  =1000 
WC 1 =-  OR  0 CFM  VARIES 
WC2=-  OR  0 ZN  VARIES 

JT=  FIRST  ESTIMATE  OF  CONCRETE  TIME  . • 

DTT=  EARTH  TIME 

DIMENSION  XWS( 150) ,XHA( 150) ,XHS(  150) ,XHW(  150 ) , XF S ( 1 50 ) , XVA ( 1 50 ) , XV 1 1 851M04 
1S( 150) ,QS( 50) » CFM ( 300 ) ,ZN(300) , 

2DBVOOO)  ,DPV(  300)  ,WBV(  300)  ,DBC(  3,300  ) ,W8C(  2,300 ) , PS  ( 10,300), 

3TVI300) ,DPC(300) , LA (300, 10) , LB( 300 , 10 ) , XDB( 300) ,EDB(100) ,EWB( 100), 

4ET( 50,50 ) , TWC ( 6 ) ,EC(6),EG(6),TS(6),TSS(6) 

5 , TGI  10,6) ,TC( 5,6) ,S (6) ,D(6) ,DD(6),ZLC6) ,DG(6) ,AG(6) , ACC  6) 11851004 

6»CKC6),CG(6) ,H(6) , HD { 6 ) , OWS ( 6 ) , QWL ( 6 ) , TQWS ( 6 ) , TQWL ( 6 ) , TQW { 6 ) 1 1 851005 

7,  DXC6)  ,CSC6) ,QW(6) ,FLX(6) ,EC6)  , AS 1 6 ) , ES ( 6 ) , TCO ( 6 ) , T I ME ( 300 ) 

COMMON  XWS*XHA,XHS»XHW»XFS»KWS»KHA , KFS * KHW, KHS»EDB*EWB,ET , KET , N11851M09 

lET,QS,C01,C02,U0tTW,GStGL,ZLEWStCFMf DPS.DBS* WBS,ETSltETS2,WCC,C03, 11851M10 
2WCD,DPA ,ZN,QG,QGL,DBV,DPV,DBC,DPC, WBC,  S , TA , DT , KWC , NN , XVA , 
3XVS,KVA,KVS,PS»LA,LB,XDB*WBV,WClfWC2  , DTT 


77 

FORMAT ( 1H010F 10.1) 

1 

FORMAT ( 10F7.0) 

T1 85 100 9 

2 

FORMAT { 1017) 

11851010 

12 

FORMAT { 43H 1 

UNDERGROUND  FALLOUT 

SHELTER  DATA 

) 

11851057 

13 

FORMAT ( 43H0 

I N IT  I 

AL  EARTH 

TEMPERATURES 

) 

11851058 

14 

FORMAT ( 43H0 

INITI 

AL  CONCRETE  WALL 

TEMPERATURE 

) 

11851059 

15 

FORMAT ( 43H 1 

TIME 

VARIABLES 

) 

11851060 

16 

FORMAT ( 43H0 

DB  ( N ) 

TV(N) 

DPV  ( N 

) QSUN(N) 

) 

11851061 

17 

FORMAT ( 4F1 0® 1 ) 

11851062 

18 

FORMATt  63H 

TG  ( 1 ) 

TG  ( 2 ) 

TG(  3 

) TG( 4 ) 

TG  ( 5 ) 

TGI 1851063 

1(6) 

) 11851064 

( 19 

FORMAT! 63H 

ICC  1) 

YC  ( 2 ) 

TCO 

) TC ( 4 ) 

TC(5) 

TC11851065 

1(6) 

) 11851066 

20 

FORMAT ( 6F10. 1 ) 

11851067 

124 

FORMAT ( 43H1 

PHYSI 

CAL  DATA  USED  FOR 

COMPUTATION 

) 

11851081 

125 

FORMAT ( 110H0S(K) 

DCK  ) 

ZL(K) 

CKCK) 

CG(K) 

ACCK11851082 

C-31 


H ( K ) 


( K ) H ( K ) HD ( K ) 

10F10.4) 

16H0  QWS(M) 

1 6H0  QWL(M) 

1 6H0  TWCtM) 

16H0  FLX(M) 

16H0  TS(M) 

14H  CONC.TEMP.  Ill, 

14H  EARTH  TEMP  111, 

120H-  EXPOSURES 
4 

100H0 
EG  ( 4 ) 


I ) AG 

126  FORMATt 

64  FORMATt 

65  FORMAT ( 

6 FORMAT ( 

67  FORMAT ( 

68  FORMAT ( 

69  FORMATt 

74  FORMATt 
63  FORMATt 

1 3 

203  FORMATt 
1 

205  FORMATt 100H0 

1 EC ( 4 ) 

204  FORMATt 10H 

210  FORMATt 50H0 

75  FORMATt 10F10.2) 

76  FORMATt 60H0 

2 TS 

7100  FORMATt 10H- 
" 62  FOR  M AT  t 1 OF 10 . 0 ) 

211  FORMATt 3F15. 1 ) 

158  FORMATt  50H-  OBSERVED 

159  FORMATt 50H-  CALCULATED 
5000  FORMAT  t 72H1 


6F15.1 ) 
6F15.1) 
6F15.1) 
6F15.1) 
-6F15.1) 
6F15.1 ) 
6F15.1) 


) 11851083 
11851084 
11851206 
11851207 
11851208 
11851209 
11851210 
11851211 
11851212 

) 


EG  t 1 ) 

EG  t 2 ) 

EGt  3 

) 

EG  ( 5 ) 

) 

ECt  1 ) 

ECt  2) 

EC  ( 3 ) 

EC  t 5 ) 

EC 

1 6 ) 

) 

6F  15« 

.3) 

DBV 

DP  V 

WBV 

CFM 

) 11851M50 

DBS 

DPS 

WBS 

RHS 

ET 

SHELTER 

SHELTER 


CONDITION 

CONDITION 


11851202 

) 11851M53 
) 11851M55 


61 


7001 

8999 


FORMATt 100 HI  TIME  QVS  QVL 

! QWST  QWLT  TWCT  QAS  QAL 

FORMATt 50H1  ADIABATIC  SHELTER  RESULTS 
FORMATt 90H0  TIME  CFM  ZN 

.DBS  WBS  EFT 

FORMATt 10F10. 1 ) 

FORMATt 1 OF  12. 2) 
jlEAD  5000 
PRINT  5000 


QGS 


DBV 


QGL 


WBV 


READ  2 , NN 


READ  2, I I , JJ,KK,KWC 

11851M33 

99 

READ  1,WC  » A I , WC 1 , WC2 

PRINT  158 

IF tWC)  3,3,4 

11851013 

3 

READ  1 , TGO 

11851M35 

DO  5 M = 1 , 5 

11851015 

DO  5 1=1,10 

11851016 

5 

TGt I ,M) =TGO 

11851017 

GO  TO  7 

1 1851M38 

4 

DO  1000  M=1 , 5 

11851M36 

1000 

READ1 , (TGt I, M), 1=1, 10) 

11 851 M3 7 

7 

READ1 » ( TCO ( M ) ,M=1,6) 

1 1851M39 

DQ  8 L=1 , 5 

1 1851M40 

DO  8 M=T,6  ~T~ 

8 

TC(L,M)=TCO(M) 

11851M41 

11 

READ  1, tS(K) ,K=1,6) 

11851043 

READ  1, (DD(K),K=1,6) 

11851044 

READ  1, (DGtKJ ,K=1,6) 

11851045 

READ  1, (Dt  K) ,K=1 ,6) 

11851046 

READ  1 , rZLTR F, K=T  ,6)  ' 

11851047 

READ  1, tAG(K) ,K=1,6) 

11851048 

READ  1, ( AC ( K ) , K=1 ,6 ) 

11851049 

C-32 


READ  1, <CK(K) ,K=i,6) 

11851050 

READ  1,  <CG(K) , K= 1 » 6 ) 

11851051 

READ  1 » ( H ( K ) ,K=1 ,6) 

11851052 

READ  1, (HD(KJ ,K=1 ,6) 

11851053 

PB=29.97 

READ  1, DT,DTT,GS,GL 

CO  1= 1.08 

C02=4.5 

C03=1050. 

U0=0. 

TW=70. 

ZLE  WS=1 . 
WCC  =- 1 . 0 
WCD=- 1.0 
CALL  PREP 


TIME (1) =6. 

DO  9 N=  1 * NN 
T V ( N ) =DBV ( N ) 

1 l 851M1 5 

9 

TIME (N) =TI ME ( N-l ) + DTT 
IF ( A I ) 5005,5005,5006 

11851M31 

5006 

READ  1,  GA,ER,CF,DATT,TW1 

5005 

CONTINUE 

PRINT  12 

11851068 

PRINT  13 

11851069 

PRINT  18 

11851070 

DO  21  1=1,10 

11851071 

21 

PRINT  20,  (TG( I » M ) ,M=i,5) 

PRINT  14 

11851073 

PRINT  19 

11851074 

DO  22  L = 1 » 5 

11851075 

PRINT  20,  ITC(L,M),M=1»6) 

11851076 

PRINT  15 

11851077 

— 

PRINT  16 
DO  5015  N= 1 » NN 

11851078 

5015 

PRINT  20 ,XDB (N) ,TV(N) ,DPV(Ni 

PRINT124 

11851085 

PRINT125 

11851086 

D0127  K = 1 » 6 

11851087 

127 

PRI NT12  6 , S ( K ) ,D(K),ZL(K),CK(K),CG(K),ACIK),AG(K),H(K), HD  I K ) 

26 

DO  28  K= 1 , 6 

11851097 

— 

DX(K)=DD(K) 
AS ( K ) =AC ( K ) 

11851100 

28 

CS(K)=CK(K) 

11851101 

OST=DT 

11851102 

DO  801  K= 1 , 6 

so  r 

E (K ) »AC ( K ) *DT/ (DOCK) *DD IK)) 
00  802  K=1 , 5 

11851176 

802 

EG(K)=AG(K)*DTT/(DG(K)*DG(K) ) 

11851184 

ECM AX=E ( 1 ) 
JK=  1 

11851M16 

DO  58  K=2 , 6 

1 1 8 5 1M 1 7 

IF ( E ( K ) - EC  M A X ) 58,58,59 

1 1 851M1 8 

59 

ECM AX=E ( K ) 

JK=K 

1 1 851M20 

58 

CONTINUE 

11851M21 

IFtECMAX-0.5 ) 91,91,92 

1 1851M22 

*2 

DST=0.4*DD( JK)*DD( JK)/AC( JK) 

11851M23 

DO  93  K= 1,6 

1 1851M24 

93 

E<K)=AC (K)*DST/(DO(K?*DD(K) ) 

91 

CONTINUE 

U851M26 
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PRINT  203 

PRINT  204, { EG ( K ) ,K=1 ,5) 

PRINT  205 

PRINT  204, ( E ( K ) , K= 1 , 6 ) 

TM  = 0 • 

DO  202  K = I , 6 

201  TSS(K)=TC(2»K) 

202  TS(K)=TC(1*K) 

DO  800  N=1,NN 

TM=TM+DTT  11851105 

X1=DBV( N ) 

X2=  WBV ( N) 

G=CFM(N) 

ZP=DPV ( N ) 

PSV=PV( ZP , PB ) 

WV=0.622*PSV/ (PB-PSV)  11851107 


CALL  XTBLU(TA1 ,5. ,60,9,QS,QGP > 

QGS= ( QGP+GS ) *ZN( N ) 

IF(KTEST) 6995, 6995,6994 

6994  CONTINUE 

PRINT  8999 , QGP , QGS 

6995  CONTINUE 

CALL  DBDPWB(TA1,DPA,WBA,WT) 

IFi A I ) 5007 , 5007 , 5008 

5008  CALL  AIRC0N(GA,G,TA1,WBA,X1,X2,TW1,ER,CF,QAS,QAL»DTAA> 
GO  TO  5009 

5007  QAS-O. 

~QAL*0. 

5009  CONTINUE 


QWST=0.  " ' 11851114 

DO  31  K= 1 , 6 11851115 

QWS(K)  =H(K)*S(K)*(TA1-TS(K) ) " 11851116 

31  QWST=  QWST  +QWSIK)  11851117 


GO  TO  (32,36), IR  11851118 

T2  ZX=QGS-QVS-QWST-QAS 


DTM=0. 


53  DTM=DTM+DST 
I R=  1 
DTA=  1 • 0 
TAl=TA-20® 


11851108 

11851111 

11851109 


30  QVS=1.08*G* (TAl-TV(N)) 


11851112 
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35  T A=  TA1 


36  DDPA=0.2 

DPA1=DPA— 10» 
I R=  1 


34  TA2  = TA1-DT A 
ZX1=ABSF { ZX  ) 

TA=(TA1*ZX2+TA2*ZX1 )/(ZXi+ZX2) 
I R=2 
TA1=TA 
GO  TO  30 


33  TA1=TA1+DTA 
ZX2=ABSF ( ZX) 
GO  TO  30 


IF(ZX)  34,35,33 


11851120 

11851121 

11851122 

11851123 

11851124 

11851125 

11851126 

11851127 

11851128 

11851129 

11851130 

11851131 


11851133 


37  PSA=PV( DPA1,PB) 


WA=0.622*PSA/ (PB-PSA) 
QVL=4780.*G*< WA-WV) 


11851135 

11851136 


CALL  DBDPWB(TA,DPA1, WBA,WT) 

IF ( A I ) 5010,5010,5011 

5011  CALL  AI RCON ( G A , G , TA  , WB A , X 1 , X2 , T W 1 , ER , CF , QAS , QAL , DT AA ) 

GO  TO  5012 

5010  QAS=0. 

QAL=0. 

5012  CONTINUE 

Q WL  T = 0.  11851137 

TWCT=0.  11851138 

DO  42  K= 1,6  11851139 

PSW=PV(TS(K) ,PB) 

WS=0.622*PSW/ ( PB-PSW)  11851141 

ZP= WA-W  S 11851142 

I F C ZP ) 39,39,40  11851143 

39  QWL ( K ) = 0 . 11851144 

GO  TO  41  11851145 

40  QWL(K)= 1061.*HD(K)*S(K)*ZP  . 11851146 

41  QWLT=QWLT+QWL  ( K ) 11851147 

TWC (K)=QWL(K)/1061.  11851148 

42  TWCT  = TWCT  + TWC { K ) ' 11851149 

CALL  XTBLU(TA,5. , 60 , 9 , QS , QGP ) 

QGL=(400.-QGP+GL)*ZN(N) 


GO  TO  (43,47) »IR  11851151 

43  ZY=QGL-QVL~QWLT-QAL 

IF(ZY)  45,  46,44  . 11851153 

44  DPA 1=DP  A 1 + DDPA  11851154 

ZY2  = ABS  F ( Z Y ) 11851155 

GO  TO  37  11851156 

45  DPA2=DPA1— DDPA  11851157 

ZY1=ABSF ( Z Y ) 11851158 

0PA-(DPAl*ZY2+QPA2*ZYl)/(  ZY 1+ZY2  ) 11851159 

I R=2 ~~  ‘ 11851160 

DP A 1 = DP  A 11851161 

GO  TO  37  \ 11851162 
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11851163 


4-6  DPA=DPA  1 
47  POA=PV{ TA, P8 ) 

RHA=100.*PSA/POA  11851165 

CALL  DBDPWB(TA,DPA, WBA,WT) 

DB=T  A 
WB= WBA 

CALL  6FT1 (DB, WB.EFSl  ) 


DO  50  K= 1 » 6 11851166 

FLX<K)= (QWS(K)+QWL(K) )/S(K)  11851167 

ES(K)=AS(K)*DST/(DX(K)*DX (K) ) 

TS(K)=2 ,*ES (K)*(FLXIK)*DX(K)/CS(K)+TSS(K)+TS(K)*(0.5*1./ES(K)-1. ) ) 

49  TC ( 1 , K) =TS ( K ) 11851172 

50  CONTINUE  11851173 

IF ( KTEST)  6985,6985,6984 

6984  CONTINUE 

PRINT  8999, ( TS ( K ) ,K=1,6) 

6985  CONTINUE 

51  DO  52  L =2 , 4 
DO  52  K = 1 , 6 

52  TC(L,K) =E(K)*(TC(L-l,K)+TC(L+l»K)+( 1 . / E ( K )-2 . ) *TC ( L , K ) ) 11851178 

DO  8002  K=  1 , 6 
8002  T cc ( K ) = TC ( 2 , K ) 


IF ( DTM-DTT ) 53,54,54 

54  DO  55  K=  1 , 5 

TG(1,K)=FN(CK(K) , CG ( K ) , TC ( 4, K ) ,TG ( 2 , K ) , DD ( K ) , DG ( K ) , 0 . ) 

55  TC(5,K) =TG<1 ,K) 

TG ( 1,6) = TC (4,6) 

TC ( 5 , 6) =TG ( 1 , 6) 

56  DO  57  1=2,9 

DO  57  K = 1 , 5 

57  TG ( I , K ) =EG(K)*(TG(I-l,K)+TG(I+lvK)+( l./EG(K)-2. )*TGC I,K) ) 
60  PRINT  61 

TSM=0» 

DO  8001  K= 1 , 6 
8001  TSM=TSM+TS { K ) 

TSM=TSM/6. 

PRINT  62  , TM,QVS,QVL,QGS,QGL,QWST ,QWLT ,TWCT  , QAS , QAL 

PRINT  159 
PRINT  76 

PRINT  75,TA,DPA,WBA,RHA,EFS1,TSM 
PRINT  6* 3 

PRINT  64, (QWS(K) ,K=1 ,6) 

PRINT  65, (QWL(K) ,K=1,6) 

PRINT  66,  (TWC(K) ,K=1,6) 

PRINT  67, ( FLX(K) ,K-1,6) 

PRINT  68,  ( TS  < K ) ,K=1,6) 

PRINT  7100 
/I  DO  70  L=2 » 5 

70  PRINT  69,  L , ( TC ( L ,K ) , K= 1 » 6 ) 

PRINT  7100 
72  DO  73  1=1,10 


11851179 
11851180 
11851181 
11851182 
1 1851M46 
11851M47 


11851186 

11851198 


11851M49 

11851203 

11851213 

11851214 

11851215 

11851216 

11851217 

11851219 

11851220 

11851221 
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11851222 


73  PRINT  74,I,(TG(I,K),K=1,5) 

PRINT  210 

PRINT  7 5 » T V ( N ) ,DPV(N) , WBV ( N ) , CFM ( N ) 

P"0  CONTINUE  11851223 

PRINT  6000 
PRINT  7000 
00  8000  N=1,NN 
X1=DBV ( N ) 

X2=WBV ( N ) 

GM=CFM{ N ) 

CALL  STRC0N(X1»X2»GM) 

8000  PRINT  7001,  T I ME ( N ) ,CFM( N ) , ZN < N ) , XI , X2 , DBS , WB S , ETS1 

CALL  SYSTEM  11851224 

END 
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SUBROUTINE  PREP 

DIMENSION  XWS( 150) ,XHA( 150) ,XHS< 150) , XHW ( 1 50 ) , XFS ( 1 50 ) , XVA ( 1 50 ) , XV 
IS t 150) , QS( 50) fCFM(300) , Z N ( 300  ) » 

2DBV(300) ,DPV( 300) , DBG (3,300 ) ,DPC( 300) , PS ( 10,300) 

3 , LA (300, 10) , LB (300, 10) ,XDB(300) ,WBV(300) ,EDB(100) ,EWB( 100 

4),ET(50»50)»WBC(2»300)»S(6) 

COMMON  XWS ,XHA,XHS,XHW,XFS,KWS,KHA,KFS,KHW,KHS,EDB, EWB, ET,  KET , N 1 1 851 M09 
1ET»QS»C01»C02»U0»TW,GS,GL»ZLEWS»CFM,DPS,DBS*W8S,ETS1»ETS2»WCC ,C03, 11851M10 
2WCD»DPA»ZN»QG» QGL ,DBV,DPV»DBC,DPC»WBC»  S , TA, DT , K WC , NN , X VA , 


3XVS, KVA,KVS,PS,LA,LB,XDB, WBV, WC1 ,WC2 
26  FORMAT ( 120H0  TIME  DBV 

1 ET  QGT  QVT  QWT 

50  FORMAT ( 1017) 

52  FORMAT ( 7F1 0« 1 ) 

1 FORMATt 10F7.0) 

2 FORMAT ( 1017) 

READ  2,KWS,NWS 

READ  1 » ( XWS ( I ) » I = 1 , NWS ) 

READ  2 , KHA , NH A 

READ  1, ( XHA ( I ) ,I=1,NHA) 

READ  2,  KHS  » NHS 
READ  1,  ( XHS ( I ) ,I  = 1»NHS) 

READ  2 , KFS  »NFS 

READ  1, ( XFS ( I ) , 1 = 1, NFS) 

READ  2 , KHW  ? NHW 

READ  1, (XHW (I ) ,1=1, NHW) 

READ  2 , K VA , NV A 

READ  1 , ( XV A ( I ) , 1 = 1 , NVA ) 

READ  2 , K VS , NVS 

READ  1, ( XVS ( I ) ,1=1, NVS) 

READ  2, KQS , NQS 

READ  1, (QS(I) ,1=1, NQS) 

READ  1, ( EDB { I ) , 1=1,29) 

READ  IT  (EWB (I) , 1=1,30) 

DO  73  1=1,29 

73  READ  1 , ( ET ( I * J ) , J= 1 , 30 ) 

IF(WCl)  4,4,5 

4 READ  1, (CFM(N) ,N=1,NN) 

GO  TO  6 

5 READ  1 , CFMO 
DO  7 N= 1 , NN 

7 CFM ( N) =CFMO 

6 I F ( WC2 ) 8,8,9 

8 READ  1 , ( ZN ( N ) , N= 1 , NN ) 

GO  TO  10 

9 READ  1 , ZNO 
DO  11  N=1 , NN 

11  ZN ( N ) =Z  NO 
10  CONTINUE 
28  FORMAT ( 12F 1 0. 2) 

TIME=0. 

STVS  ( 1 ) +S (2 ) +s ( 3 ) +S ( 4 ) + S (5 ) +S (6 ) 

IF(KWC)  61,61,  62 

61  DO  300  N=1 , NN 
READ  50,(LA^N,L),L=1,10) 

R E A D 50  , ( LB ( N , L ) , L = 1 , 1 0 ) 

GO  TO  500 

62  READ  1 , ( DB V(  N) ,N= 1 , 240 ) 

READ  1 , ( WB V ( N) ,N=1*240) 

READ  1, (DBC(i,N),N=l, 240 ) 


, DTT 

WBV  DBS  WB  S 

QWS  TSM  HTX ) 

11851002 

11851024 

11851011 

11851012 

11851015 

11851016 

11851017 

11851018 

11851019 

11851020 

11851021 

11851022 

11851023 

11851024 

11851025 

11851026 

11851027 

11851028 

11851029 

11851030 


11851048 


11851018 

11851F05 


11851FQ9 
1 1851F10 
1 1851F12 
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READ  l» ( WBC ( 1 » N ) » N= 1 , 240 ) 

READ  1, (DBC(2,N) ,N=i,240) 

READ  1, (WBC(2,N) ,N=1,240) 

READ  1, (DRC(3,N) ,N=1,240) 

READ  If ( XD8 ( N ) ,N=1, 240 ) 

D063  J=l,10 

63  READ  1, (PS< J,N) ,N=1,240) 

DO  600  N= 1 ,240 
600  CFM(N)=CFM<N)/2. 

500  PRINT  26 

DO  400  N= I , NN 
IF(KWC)  401,401,402 
401  DBV(N)=LA(N,1 )/10 
X=DB  V ( N ) 

WBV(N)=LA(N»2)/10 
Y=WBV(N ) 

CALL  DB WBDP ( X , Y , Z , WQ ) 

DPV ( N ) = Z 

X=( LA(N,3) +LA(N,5)+LA(N,7)+LA(N,9) )/40 
DBSM=X 

Y=(LA(N,4)+LA(N,6)+LA(N»8)+LA(N,  10)  )/40 
WBSM=Y 

CALL  DBWBDP(X,Y,Z,WQ) 

DPSM=Z 

X=(LB(N,l)+LB{N,2)  + LB(N,3)+LB(N,4)-*-LB(N,5)+LB(N,6)+LB(N,7)+LB(N, 
1 +LB ( N » 1 0 ) )/90 
I F ( L B ( N » 4 ) ) 53,53,54 


1 1 85  IF  1 3 
1 1851F14 
11851F15 


11851F18 

11851F19 


11851014 


11851010 

11851012 

11851014 

11851016 

8)11851018 

11851019 

11851020 


53  X=X*9. / 8 . 

54  T SM=  X 

XDB(N)=LB(N,9)/10 

GO  TO  70  11851F07 

402  X=DB V ( N ) 

Y=WB V { N ) 

CALL  DBWBDP(X,Y,Z ,WD)  11851F23 

DPV ( N ) = Z 

X=  ( DBC  ( 1,N)+DBC(2,N)+DBC(3,N)  )/3. 

DBSM=X 

Y=<WBC< 1,N)+WBC(2,N) )/2. 

WBSM=Y 

CALL  DBWBDP( X, Y,Z,WC)  11851F27 

DPSM=Z 


SUM=0. 

SUN=0. 

DO  64  J=l,10 

JF(PS(J,N))  64,64,66 

66  SUM= SUM + PS ( J , N ) 
SUN=SUN+1. 

64  CONTINUE 

67  TSM=SUM/SUN 
70  CONTINUE 

DPC ( N)=DPSM 
T I ME-TI ME+DTT 
X=DB V { N ) 

Y=DP V ( N ) 

Z=WB V ( N ) 

CALL  TBLU( X,KWS,XWS,WSS) 
CALL  TBLU(Y,KWS,XWS,WSV) 
CALL  TBLU( X,KVA,XVA,VA1) 
CALL  TBLU< X , K VS , XVS , VS 1 ) 
V=VAl+( VS1-VA1)»WSV/WSS 


11851F31 
1 1851F32 
I1851F33 
11851F34 
1 1 851F35 
1 1851F36 
11851F37 

1 1851F08 


11851106 

11851107 

11851108 

11851109 

11851110 
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XX=DBSM 
YY=DPSM 
ZZ= WBSM 

CALL  EFT1 ( XX, ZZ,SET) 

CALL  TBLU( YY 9 KWS , XWS , WS A ) 

CALL  XTBLU (XX, 5. ,60,9,QS, QGS ) 

QGT=(40Q.+GS+GL)*ZN(N) 

QGS  =QG  S*ZN ( N > 

WGL  =(400. -QGS ) *ZN { N ) / 1053 . 7 
QVS=CFM(N) /V*0. 24*60. MXX-X) 

WVL=CFM ( N ) /V*60. * < WSA-WSV ) 

QVT=QVS+WVL*1054. 

QWS=QGS-QVS-GS 
QWT=QGT-QVT 
HTX=QWT/ST 
RATI 0=Q WT/QGT 

PRINT  28»TIME»DBV(N) » WBV ( N ) , DBSM , WBSM , SET , QGT , QVT » QWT , QWS , TSM , 
1HTX 

400  CONTINUE 

T A=DBC ( 1,1) 

DP A=DPC ( 1 ) 

RETURN 

END 


11851024 

11851113 

11851114 

11851118 

11851117 

11851115 

11851116 

11851026 

11851027 

11851028 


11851029 


c 


11851E02 


SHELTER  AIR  CONDITION  CALCULATION 
SUBROUTINE  STRCON ( DB, WB* Z FM ) 

DIMENSION  XWS! 150) ,XHA{ 150) ,XHS( 150) ,XHW( 150) ,XFS(150) , 

1EDB(  100)  ,E'WB(  100)  ,E T ( 50 , 50 ) , QS  ! 50  ) ,CFM(300) 

COMMON  XWS,XHA,XHS,XHW,XFS,KWS,KHA,KFS,KHW,KHS, EDB,  EWB  » ET  » 
3KET,NET,QS,C01 , C02 , UO , TW , GS , G L, ZL EWS , CFM , DP S , DBS , WB S, ETS 1 , ETS 2 
4wCC,C03, WCD 


12 

FORMAT ( 10H 

ETS  1 

F 10 

13 

FORMAT! 10H 

ETS2 

F 10 

9 

FORMAT! 10H 

QGL 

F 10 

4 

FORMAT! 10H 

DBS 

F 10 

5 

FORMAT! 10H 

QGS 

F 10 

14 

FORMAT! 10H 

DPS 

F 1 0 

6 

FORMAT! 10H 

TWS 

F 10 

7 

FORMAT! 10H 

DP 

F 10 

11 

FORMAT! 10H 

WBS 

F 10 

10 

FORMAT! 10H 

WA 

F 10 

8 

FORMAT! 10H 

WV 

F 10 

2) 


Y 1 = C01*  ZFM+UO 
Y2=C01*ZFM*DB+U0*TW 
Y3=Y2/Y 1 

I F ( Y3-100. ) 30,30,31 
30  X 1 = Y 3 
3 X2=X 1+0 • 5 

CALL  XTBLU(X1,5. ,60,9, QS, QS1 ) 
CALL  XTBLUIX2 ,5. ,60,9,QS, QS2) 
Zl=( QS1+GS )/Yl+Y3-Xl 
Z2=!QS2+GS)/Y1+Y3-X2 
I F ( Z 2 ) 1,1,2 
? X1=X2 


GO  TO  3 
1 Z2=ABSF(Z2) 

DBS=X1+0.5*Z1/Z2/ ( 1.+Z1/Z2) 

CALL  XTBLU(DBS,5. ,60, 9 , OS , QGS ) 

GO  TO  32 

31  DBS=( 1400.+GS+Y2)/! Yl+14. ) 

QGS=-14.*( DBS-100. ) 

32  CONTINUE 

CALL  TBLU( TW,KWS,XWS,TWS) 

CALL  DBWBDP!DB,WB,DP,WV) 

QGL=AOO.-QGS 

ZF=ZLEWS*U0/0.243 

WA=( (QGL+GL ) /C03+ZF*TWS+C02*ZFM*WV)/ (C02*ZFM+ZF ) 
IF(WCC)  19,20,20 

20  PRINT  10, WA 
19  CONTINUE 

IF(WA-TWS)  101,101,102 

101  IF ( WCD)  103,103,102 

103  WA=WV+< GGL+GL) / ( C02*ZFM*C03 ) 

102  CALL  ,DBWWBH(DBS,WA,WBS,H) 

IF(WCC)  21,22,22 

22  PRINT  1 1 , WBS 

21  CONTINUE 

CALL  DBWBDPI DBS , WBS , DPS , WA ) 

IF(WCC)  23,24,24 
PRINT  14, DPS 
PRINT  1 0 , WA 

23  CONTINUE  ’ " 

IF(WCC)  25,16,16 

16  PRINT  4 , DB S 


1 1851E05 
,11851028 


11851 


11851E10 
1 1 85 IE  1 1 


1 1851E14 


11851E15 

11851E16 

11851E17 

11851 

11851 

11851 

11651 

11851E 1 8 

11851E19 

11851E20 


1I851E24 
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PRINT  5 , QGS 
PRINT  6 ,TWS 
PRINT  7,  DP 
PRINT  8,WV 
PRINT  9 « QGL 
PRINT  10, WA 
PRINT  1 1 » WBS 
PRINT  14, DPS 
25  CONTINUE 

CALL  EFT1( DBS, WBS  »ETS1 ) 
IF(WCC)  15,17,17 
17  PRINT  12 , ETS 1 
15  RETURN 
END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


2 

3 

99 


100 


FOR 

SUBROUTINE  A I RCON ( GA »RW »SP  * CF  *GW * DBS t WBS  * TW I > Q AS » QAL * NROW ) . 

GA  = A I R FLOW  THROUGH  COIL  »CFM 

GW  = COOL ANT  FLOW  ( LBS/HR ) * ( SPEC  I F I C HEAT)  * 3TU/HR/DEG=500*GPM  FOR 

R W = C 0 1 L THERMAL  RESISTANCE  EXCLUDING  AIR  FILM  RESISTANCE 

SP=CO I L AIR  SIDE  HEAT  TRANSFER  SURFACE  PER  ROW 

CF=CO I L CONTACT  FACTOR  PER  ROW 

NROW  = COIL  ROWS 

DBS  » WBS  » =C0 1 L INLET  AIR  DRY-AND  WET-BULB  TEMPERATURES 
TWI  = CO  I L COOLANT  ENTERING  TEMPERATURE 

QAS  =C0 I L SENSIBLE  CAPACITY 
QAL  =C0 I L LATENT  CAPACITY 

DIMENSION  XWS  (150)  »XHA  ( 150  ) »XhS  ( 150  ) »XHW  (-150  )»XFS(150)*T(20)» 
1 ) *TW( 20 ) *TAL ( 2 ) * TWG ( 2 ) »FGS ( 2 ) 

COMMON  XWS , XH A » XHS  » XHW » XFS  » KWS * KHA » KFS  * KHW » KHS 
T ( 1 ) =DBS 

CALL  DBWBDP ( DBS»WBS »DPS , WX ) 

W ( 1 ) = wx 

TWG ( 1 ) = D B S - 1 • 

TWG( 2 ) = T W G ( 1 )-l . 

A = SP / ( RW*GA*C  F *4 • 5 ) 

7 DO  1 00  K = ] » 2 
TW ( 1 ) =TWG ( K ) 

DO  99  I = 1 » NROW 
X = T ( I ) 

Y = W ( I ) 

Z = TW(  I ) 

CALL  SURFT( A*X»Y*Z»TSX*WSX ) 

T(  1 + 1 )=T( I )-(T ( I )-TSX)*CF 
Z=WSX-W( I ) 

I F (Z ) 1 »2  *2 

1 W ( I +1  ) =W ( I ) - ( W ( I ) -WSX ) *CF 
GO  TO  3 
W(  I+l  ) = W ( I )' 


T W ( I + l ) = T W ( I ) - 
CONTI NUE 
NZ=NROW+ 1 
FGSCK ) =TW ( NZ ) -TWI 
TAL ( K ) =T ( NZ ) 

T EST=  FGS ( 1 ) *FGS ( 2 ) 
IF(TEST)  4»5»6 
T W G ( 1 )=TWG( 2 ) 

T WG ( 2 ) =TWG ( 2 ) ~ 1 • 

GO  TO  7 


( SP / (GW*RW)  ) * ( T SX-TW ( I ) ) 


11851001 
11851002 
11851003 
H2011851004 
11851005 
11851006 
11851007 
11851008 
11851009 
11851010 
11851011 
11851012 
W ( 2011851013 
11851014 
11851015 
11851016 
11851017 
11851018 
11851019 
11851020 
11851021 
11851022 
11851023 
11851024 
11851025 
11851026 
11851027 
11851028 
11851029 
11851030 
11851031 
11851032 
11851033 
11851034 
11851035 
11851036 
11851037 
11851038 
11851039 
11851040 
11851041 
11851042 
11851043 
11851044 


5 TWL=T WG ( 2 ) 

TL=TAL ( 2 ) 

GO  TO  98 

4 C = ABS F ( FGS ( 1 ) / FGS ( 2 ) ) 

TWL= ( TWG ( 1 ) + T W G ( 2 ) *C ) / ( 1 . + C ) 
TL= ( TAL( 1 )+TAL ( 2 ) ) / 2 . 

98  QAS=1 .08*GA*( DBS-TL ) 

Q A L = G W * ( T W L — T W I ) —QAS 

RETURN 

END 


11851045 
11851046 
11851047 
11851048 
11851049 
11851050 
11851051 
1 185.10  52 
11851053 
11851054 


C -43 


c 


FOR 

SUBROUTINE  SURFT( A»TA»WA*TW*TS*WS) 

SOLVE  FOR  TS  FROM  0. 243* ( TA-TS ) +1060 .*( WA-WS ) =A* ( TS-TW ) 
DIMENSION  XWS ( 150 ) >XHA ( 150  ) ,XhS ( 1 50 ) *XHW ( 1 50 ) > XFS ( 150 ) * Y ( 2 ) » F ( 2 ) 
COMMON  XWS»XHA  »XHS  »XHW  »XFS  *KWS*KHA*KFS »KHW*KHS 
Y ( 1 ) =T A 
Y ( 2 ) =TA-1. 

5 DO  1 1=1 >2 

P1=PV ( Y ( I ) *29. 92  ) 

W = 0.622*P1/ ( 29 .92-PI  ) 

1 p ( I ) =A*( Y( I ) -TWJ-0.243* ( TA-Y ( I ) ) +1060.* ( WA-W ) 

T EST  = F ( 1 )*F ( 2 ) 

IF(TEST)  2>3>4 
A Y ( 1 ) = Y ( 2 ) 

Y(2)=Y(2)-1. 

GO  TO  5 
3 TS= Y ( 2 ) 

GO  TO  10 

2 C = ABSF ( F ( 1 ) / F ( 2 ) ) 

T S = ( Y ( 1 )+Y( 2 )*C ) / ( 1 .+C) 

10  P l =PV ( TS, 29.92  ) 

WS=0 . 622*P1 / ( 29.92-Pl ) 

return 

END 


11851001 

11851002 

11851003 

11851004 

11851005 

11851006 

11851007 

11851008 

11851009 

11851010 

11851011 

11851012 

11851013 

11851014 

11851015 

11851016 

11851017 

11851018 

11851019 

11851020 

11851021 

11851022 

11851023 

11851024 
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FUNCTION  FN(XI,X2tYI,Y2,Zl*Z2.Z3)  11853720 
Pl=  X1*Y1/Zl  +X2*Y2/Z2  + Z3  11853730 
P2=  XL/Z1  + X2/Z2  1 1853740 
**N  = P 1 /P  2 11853750 
RETURN  11853760 
END 


C-45 


EFFECTIVE  TEMPERATURE  CALCULATION 
SUBROUTINE  E FT  1 ( DB , WB , EF X ) 

DIMENSICN  XWS ( 150 ) , XHA ( 1 50 ) f XHS ( 150) ,XHW( 150 ) ,XFS< 150  ) , 

EDB ( 100 ) ,EWB ( 100 ) ,ET<  50,50) 

. .OMMON  XWS , XHA, XHS, XHW, XFS.KWS, KH A, KFS,KHW,KHS, EDB, EWB,ET, KET.NET 
IF(DB-WB)  1,1,2 

1 EFX=DB 
GO  TO  15 

2 I F ( DB-ECB ( 29 ) ) 3,4,4 

3 EFX=0. 


GO  TO  15 

4 I F ( WB-E WB ( 30  ) ) 3,5,5 

5 I F ( DB-ECB ( 1 ) ) 6,6,7 

7 EFX  = 10C. 

GO  TO  15 

6 IX=0 


DO  8 1=1,29 

I F ( DB-ECB ( I ) ) 9 

10  X 1=EDB ( I ) 

X2  = EDB( 1-1  ) 

GO  TO  11 

9 IX=IX+1 
8 CONTINUE 

11  J Y=0 

DO  12  J = 1 » 30 
I F ( WB-E  WB ( J ) ) 13,14,14 

14  Y l = E WB ( J ) 

Y2=EWB( J-l  ) 

GO  TO  16 


1®iY=JY+1 

1?»C0NTINUE 

16  Z 1 1=ET{ IX+1, JY+1) 


,10,10 


Z 12  = ET { IX, JY+1  ) 

Z22  = ET(  IX, JY) 

Z2 1=ET ( IX+i, JY) 

IF [111* 112*122*121 ) 3,3,17 
17  Z2=Z22+(DB-X2)*(Z21-Z22)/ (X1-X2) 
Z1=Z12+ (DB-X2 )*(Z11-Z12)/ (X1-X2) 
EFX  =Z2+(WB-Y2)»( Zl-Z2)/( Y1-Y2) 
15  RETURN 
END 


11851  44 


11851  45 


C-46 


enthalpy  calculation  by  cb  and  wb 

1 1851G02 

SUBROUTINE 

DBWBH(DBfWBtH) 

11851G03 

DIMENSION  XttSI 150) , XHA( 150) »XHS( 150) « XHWt 150  ) , XFS ( 150 ) 

1 1851S04 

COMMON 

1 XWS 

, XHA ,XHSf  XHW,XFS , KWS , KHA , KFS , KHW,KHS 

11351G05 

;all 

TBLU 

(WB»KWS»XWS»WS) 

1 1 851G06 

CALL 

TBLU 

(WBtKHWfXHWtHW) 

11851G07 

CALL 

TBLU 

(DB,KWS,XWS,WSS) 

11851G08 

CALL 

TBLU 

( DB,KHA , XHA, HA  ) 

1 1 851G09 

CALL 

TBLU 

( WB,KHS ,XHS,HS ) 

11851G10 

CALL 

TBLU 

(DB,KHS,XHS,HSS) 

11851G11 

HAS=HSS-HA 

1 1 8 5 1G  l 2 

X=(HS- 

HI**WS)*HAS-HA*HW*WSS 

11851G13 

Y=HAS- 

Hfe*WSS 

11851G14 

H = X / Y 

11851G15 

RETURN 

1 1 851G16 

END 


C-47 


RELATIVE  HUMIDITY  CALCULATION  BY  DB  AND  WB 
SUBROUTINE  DBWBRH ( DB , WB , PB, RH ) 

DIMENSION  XWS( 150) ,XHA( 150),  XHS ( 150 ) ,XHW ( 1 50 ) , XFS ( 150 ) 
COMMON  XWS,XHA, XHS, XHW,XFS , KWS, KHA, KFS,KHW,KHS 
^ ALL  DBWBOP(DBfkB»DPfW) 

CALL  TBLU(CB»KtoS»XWS,WS) 

ZM=W/WS 
DBB  = DB/  10. 

CALL  TBLUl DBB,KFS»XFS»FS) 

Z=PV(DB,P8)/PB 

RH= 100. *ZM/ { 1 . — ( 1 . -ZM ) *FS  ftZ ) 

RETURN 

END 


11851B1 L 
11851B12 

11851B14 

11851B17 

11851B16 


1 1851B20 
1 1851B2  1 
11851B22 
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WET  BULB  CALCULATION  BY  DB  AND  W 11851F02 

SUBROUTINE  DBW WBH ( DB , W , WB , H ) 

DIMENSION  XWS< 150) ,XHA( 15  0) ,XHS(  150) ,XHW(  150) ,XFS( 150) 

COMMON  XWS,XHA,XHS,XHW,XFS,KWS,KHA,KFStKHW,KHS  11851F05 

CALL  TBLU(DBfKHA,XHAfHA)  11851F06 

CALL  TBLU( DB,KHS ,XHS,HS  ) 11851F07 

CALL  TBLU(DB,KWSfXWS,WS) 

IF(W-WS)  4,5,5 
5 HB=DB 
H=HS 


GO  TO  6 

4 H=HA-»-W*(HS-HA)/WS 


X 1 = DB 

11851F10 

3 

X2=X 1-0 . 5 

1 1 85  IF  1 1 

CALL  TBLU(X1,KHS,XHS,HS1) 

11851F12 

CALL  TBLU(X2,KHS,XHS,HS2) 

1 1851F13 

CALL  TBLU(X1,KWS,XWS,WS1) 

11851F14 

CALL  TBLU(X2,KWS,XWS,WS2) 

11851F15 

CALL  TBLU(X1,KHW,XHW,HW1) 

1 1851F16 

CALL  TBLU( X2 , KH W , XHW , HW 2 ) 

1 1 85  IF  1 7 

Z1=HS1-HW1*(WS1-W) 

11851F18 

Z2=HS2-HW2*(WS2-W) 
I F ( Z2— H ) 1,1,2 

11851F19 

2 

X 1 = X 1—0  • 5 

11851F21 

GO  TO  3 

11851F22 

1 

WB=X2+0.5* (H-Z2 ) / ( Z1-Z2 ) 

6 

RETURN 

END 

C-49 


VAPOR  PRESSURE  CALCULATION 
FUNCTION  P V ( X » PB ) 

T=X+459.688 
T S=671 .688 

fMl=-( 7.90298 }*( TS/T-1. ) 

TM2=(5.C2808)*(L0G10F(TS/T) ) 

TM3=-( I .3816/ ( 10.** <7. )>)*(< 10.** ( I I . 344* ( 1 . -T/TS )))-!.) 
TM4=  (8.  132  8/1000.  )*<  10.  **  (-3.49149*  ( TS'/T-l.  ) >-l  . ) 

ANS1=TM 1+TM2+TM3*TM4 
ANS2=10.»* ( ANSI ) 

PV=PB*ANS2 

RETURN 

END 


11851H1* 
1 1 8 5 1H l ‘ 
11851HH 
1 1851H1', 
11851HH 
1 1 85  1H 1 c 
11851H2C 
11851H21 
11851FI2; 
L1851H22 
1 1 85 1H2^ 
1 1 8 5 1H2  f 
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DEW-POINT  CALCULATION  BY  DB  AND  WB 
SUBROUTINE  D BWEDP ( DB , WB , DP , W ) 

DIMENSION  XWS( 150 ) ,XHA( 150 ) , XHS ( 150 ) ,XHW( 150 ) ,XFS( 150 ) 
COMMON  XWS,XHA,XHS,XHW,XFS,KWS,KHA,KFS,KHW,KHS  ' 
.F(DB-WB)  1,1,2 

1 DP=DB 

CALL  TBLU(DB,KWS,XWS»W) 

GO  TO  100 

2 CALL  TBLU(CB,KHA,XHAVHA) 


5 IS  FULL.  GIVE  ME  ANOTHER  A3  AND  START. 

PERMAMENT  TAPE  REDUNDANCY  WAS  DETECTED  WHILE  READING  THE  RECORD  THAT  PRODUCED 


CALL  TBLU(DB,KFS,XHS»HS) 
CALL  TBLU (CB»KWS,XWS,WS) 

1185101} 

1185101; 

CALL  TBLU (WB  ,KHS  , XHS, HSWB) 

CALL  TBLU( WB, KHW  ,XHW,HW) 

CALL  TBLU ( WB,KWS  ,XWS,WSWB) 
W=WS*(HSWB-HA-HW*WSWB) / ( H S-H A-HW*WS ) 

1185101: 
1185101" 
1 1 8 5 1 0 1 f 
118  5101 < 

5 

Y 1 = WB 
Y2= Yl-0 • 5 

1185101' 
1185101 J 

CALL  TBLU( Y1 ,KWS,XWS  ,Z1 ) 
CALL  TBLU( Y2,KWS,XWS,Z2) 

1185101* 
1 1851021 

I F ( Z2— W ) 3 » 3 , 4 

1185102 

4 

Y 1 = Y 1-0  • 5 

1185102; 

GO  TO  5 

1185102 

3 

DP=Y2+0.5* (W-Z2)/ (Z1-Z2) 

1185102- 

100 

RETURN 

END 

1185102 

11851CC; 
11851C0: 
1185100" 
1 1851005 
1185100* 
1185100' 
1185100J 
1185100" 
118510K 
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toE  T BULB  CALCULATION  BY  DB  AND  DP  1I851BO 

SUBROUTINE  DBDPWB ( DB , DP , WB , W ) LI851BO 

DIMENSICN  XWS ( 150 > ,XHA( 150) , XHS ( 150 ) , XHW ( 1 50 ) , XF S I 1 50 ) 

COMMON  XWS,XHA,XHS,XHW,XFS,KWS,KHA,KFS,KHW,KHS  11851B0 

CALL  TBLU(DP,KWS»XWS,W)  11851B0 

CALL  DBfcWBHI DB ,WvViB»H) 

RETURN  11851B0 

END 


C-52 


LINIER  INTERPOLATION  1I85IH0 

SUBROUTINE  TBLU ( X , KX , V ♦ Y ) 1I851H0 

DIMENSION  V ( 150) 

I = X I N T F ( X ) 1I851H0 

4_= J-KX+  2 11851H0 

LL= J-KX+ 1 11851H0 

VU=V ( L ) I185IH0 

VL= V ( LL ) 11851H0 

Y= ( VU-V  L)*(X— INTF(X) )+VL  11851H1 

RETURN  1L85IH1 

END 


C-53 


SUBROUTINE  XTBLU(X»DX,KX,LX,V,Y) 
DIMENSION  V( 50) 

J=XINTF (X) 

M=KX+XINTF(DX)*(LX-1) 

IF(J-KX)  1,2,2 

1 Y=V(1) 

GO  TO  1C 

2 I F ( J— M ) 3,4,4 

4 Y=V ( LX ) 

GO  TO  1C 

3 DO  6 1=1, LX 
Z=I 

P = K X 

Q=P+Z*DX 

IF(X-Q)  5,5,6 

5 Q 1 = V ( 1+1 ) 

Q2=V ( I ) 

Y=Q1-(Q1-Q2)*(Q-X)/DX 
GO  TO  1C 

6 CONTINUE 
10  RETURN 

END 


11851C0 

1 1 85 ICO 
1 1 85  ICO 

1 1 85  ICO 
1 1 85  ICO 
1 1 8 5 ICO 
11851C1 
11851C1 
1 1 85 1C  1 
11851C1 
1 1851C1 
1 1 85  1C  1 
1 1 8 5 1C  1 
1 1 8 5 1C  1 
11851C1 
1 1 8 5 1C  1 

11851C2 

11851C2 


C-54 


Sample  Input  and  Output  Data  on  M-4  Program 


UNDERGROUND  FALLOUT  SHELTER  DATA 


INITIAL 

c ART  H 

TEMPERATURES 

T G ( L ) 

T r.  ( 20 

TG  ( 3 ) 

TGI  4) 

TGI  5) 

TG  I 6 ) 

70.5 

72.  0 

72. C 

72.0 

73.0 

70.0 

72.  0 

72. C 

72.0 

72.0 

69.0 

71.  0 

71.  C 

71.0 

72.0 

67.0 

70.  0 

70. C 

70.0 

71.0 

67.0 

70.  0 

7C.0 

70.0 

71.0 

67.0 

70.  0 

7C.0 

70.0 

70.0 

67.0 

70.  0 

7C.C 

70.0 

69.0 

67.0 

69.  0 

69. C 

69.0 

68.0 

67.0 

68.  0 

68.  C 

68.0 

67.0 

67.0 

67.  0 

67. C 

67.0 

67.0 

. ..  ... 

INITIAL 

r0NC  R E T E WALL 

TEMPERATURE 

TCI  1 ) 

T r ( 2 ) 

TC  ( 3 ) 

TC  I 4 ) 

TCI  5) 

TCI6) 

73.0 

73.  0 

73.0 

73.0 

73.0 

73. 

73.0 

73.  0 

73. C 

73.0 

73.0 

73. 

73.0 

73.0 

73. C 

73.0 

73.0 

73. 

73.0 

73.  0 

73. C 

73.0 

73.0 

73. 

73.0 

73.  0 

73. C 

73.0 

73.0 

73. 

0-1 


f I ME  VARIABLES 


OB  ( N ) 

T'M  N) 

70.0 

76. 

72.0 

72. 

74.0 

74. 

75.0 

66 . 

77.0 

68. 

77.0 

66. 

76 . 0 

64. 

76.0 

63. 

75.0 

62. 

75.0 

81. 

75.0 

91. 

75.0 

92. 

77.0 

93. 

7 8.0 

92. 

79.0 

88. 

80.0 

85. 

80.0 

81. 

80.0 

79. 

80.0 

77. 

80.0 

76. 

80.0 

79. 

81.0 

84. 

81.0 

87. 

80.0 

90. 

82.0 

92. 

82.0 

92. 

83.0 

88. 

83.0 

84. 

83.0 

82. 

83 . 0 

80. 

8 3.0 

79. 

82.0 

78. 

82.0 

80. 

82.0 

86. 

82.0 

91. 

82.0 

92. 

03.0 

92. 

83.0 

92. 

84.0 

88. 

84.0 

83. 

84.0 

81. 

84.0 

79. 

84.0 

80. 

84.0 

79. 

84.0 

81. 

83.0 

86. 

83.0 

90. 

83.0 

92. 

84.0 

92. 

84.0 

92. 

85.0 

88. 

85.0 

8 3. 

85.0 

80. 

85.0 

80. 

84.0 

79. 

84.0 

79. 

84.0 

79. 

CPV(N)  QSUN(N) 

57.2 
47  . 9 
61.9 

48.6 
55. C 

48.6 

48.2 
47.  C 

43.5 

68 . 2 

71.8 
71.  A 

71 . C 

68.2 

69.9 

71. 1 

71.3 

69 . 1 

68.  A 

71.8 
7C.6 
73  . C 

71.9 

72.2 

71.4 
73. C 
73. C 

71.6 

72.3 
7 C . 2 

69.  1 

69.5 
7C.2 

69.2 

71.8 

71.4 

69.8 

71.4 

71.4 

70.4 

69.8 
69  . 1 

71.7 
7C  • 6 

71.3 

72.2 
7C.6 

69.8 

71.4 

71.4 

71.4 
7 C . 4 

70.2 
7C.2 
72.1 

72.  1 
72.  1 

0-2 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


/ 


83.0  87.0  71.9 

83.0  90.0  72.2 

83.0  92.0  69.8 

84.0  92.0  71. A 

84.0  92.0  71.4 

85.0  88.0  73. C 

85.0  83.0  72. C 

85.0  82.0  7C.8 

85.0  80.0  7C.2 

85.0  79.0  7C.6 

85.0  79.0  70.6 

85.0  80.0  73.1 

85.0  86.0  72.2 

84.0  89.0  72.6 

84.0  92.0  71.4 

85.0  93.0  71. C 

85.0  91.0  70.2 

85.0  87.0  70.3 

85.0  83.0  70.4 

85.0  81.0  69.8 

85.0  80.0  70.2 

85.0  79.0  70. 6 

85.0  78.0  71. C 

85.0  80.0  71.7 

84.0  85.0  72.6 

84.0  89.0  72.6 

83.0  91.0  71.8 

84.0  93.0  71. C 

83.0  92.0  71.4 

85.0  88.0  71.4 

86.0  83.0  70.4 

86.0  81.0  71.3 

86.0  79.0  70. 6 

85.0  79.0  70.6 

85.0  79.0  70. 6 

85.0  80.0  71.7 

85.0  82.0  73.8 

85.0  90.0  72.2 

85.0  92.0  71.4 

85.0  92.0  71.4 

85.0  91.0  71.8 

85.0  87.0  70.3 

85.0  83.  0 72- G 

85.0  80.0  71.7 

85.0  79.0  72.1 

85.0  78.0  71. C 

84.0  77.0  71.4 

84.0  79.0  72.1 

84.0  85.0  72.6 

84.0  39.0  72.6 

83.0  91.0  73.4 

84.0  92.0  71.4 

85.0  91.0  71.8 

85.0  88.0  69.9 

86.0  83. 0 72.0 

85.0  80.0  71.7 

85.0  79.0  7C.6 

85.0  78.0  71.0 

85.0  78. 0 71. C 

84.0  79.0  73.5 

0-3 


84.0 

85.  0 

72.6 

84.0 

89.  0 

72.6 

84.0 

92.  0 

73. C 

84.0 

92.  0 

71.4 

85.0 

92.  0 

73.0 

85.0 

88.  0 

73.0 

85.0 

84.  0 

7C.C 

85.0 

81*0 

71.3 

85.0 

80.  0 

70.2 

85.0 

78.  0 

71. C 

84.0 

78.  0 

71. C 

85.0 

80.  0 

71.7 

84.0 

85.  0 

72.6 

84.0 

90.  0 

72.2 

84.0 

92.  0 

71.4 

84.0 

92.  0 

71.4 

85.0 

91.  0 

71.8 

85.0 

88.  0 

69.9 

86.0 

83.  0 

72. C 

86.0 

81.0 

71.3 

85.0 

79.  0 

70.6 

86.0 

79.  0 

7C.6 

85.0 

78.  0 

71. C 

85.0 

79.  0 

72.1 

86.0 

85.0 

72.6 

85.0 

89.  O' 

72.6 

85.0 

90.  0 

72.2 

85.0 

92.  0 

73.0 

86.0 

91.0 

71.8 

86.0 

88.  0 

69.9 

86.0 

84.  0 

7C.C 

86.0 

81.  0 

71.3 

86.0 

78.  0 

71.0 

86.0 

78.  0 

72.5 

86.0 

78.  0 

71. C 

85.0 

80.  0 

71.7 

85.0 

85.  0 

72.6 

85.0 

88.  0 

73.0 

84.0 

92.  0 

71.4 

85.0 

93.  0 

72.6 

85.0 

92.  0 

71.4 

86.0 

88.  0 

69.9 

86.0 

84.  0 

71.6 

86.0 

81.  0 

71.3 

86.0 

79.  0 

69. 1 

86.0 

79.  0 

72.  1 

85.0 

79.  0 

72.1 

85.0 

80.  0 

71.7 

85.0 

85*  0 

71.1 

85.0 

90.  0 

72.2 

85.0 

92.  0 

71.4 

85.0 

92.  0 

71.4 

86.0 

'92.  0 

71.4 

86.0 

89.  0 

71.0 

86.0 

83.  0 

72. C 

86.0  80.0  71.7 

86 . 0 78.7) 71.0 


86.0 

78.  0 

71.0 

86.0 

77.  0 

71.4 

85.0 

80.  0 

73.  1 

0-4 


85.0 

85.0 

85.0 

86.0 

86.0 

86.0 

85.0 
6 5.0 

84.0 

84.0 

84.0 

33.0 

83.0 

82.0 

83.0 

83.0 

63.0 

83.0 

63.0 

83.0 

83.0 

83.0 

83.0 

83.0 

82.0 

82.0 


36.  0 
90.  0 

92.  0 

93.  0 
93.  0 
88.  0 

84.  0 

81.0 

79.  0 
79.  0 

79.  0 

80.  0 

85.  0 
90.  0 
92.  0 
92.  0 

91.0 
88.  0 
83.  0 

81.  0 

80.0 
79.  0 
79.  0 

81.0 

86.  0 

90.  0 


72.2 

72.2 
71.4 
71  .C 

72.6 
71.4 

71.6 

71.3 
7Q.6 

6 9.1 
7G.6 

71.7 

71.1 
7C  . 6 

71.4 

71.4 

71.8 

71.4 
72. C 
71.3 

71.7 

7 C . 6 

72. 1 

72.7 

72.2 

72.2 


0-5 
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Four  digital  computer  programs  for  simulating  thermal  Four  digital  computer  programs  for  simulating  thermal 

environment  in  underground  fallout  shelters  are  environment  in  underground  fallout  shelters  are 

described.  The  calculated  temperature  and  humidity  described.  The  calculated  temperature  and  humidity 

of  seven  prototype  shelters  are  compared  with  the  of  seven  prototype  shelters  are  compared  with  the 

experimentally  observed  values.  I experimentally  observed  values. 


